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ABSTRACT 

Adopting the framework of the halo occupation distribution (HOD), we investigate the ability 
of galaxy clustering measurements to simultaneously constrain cosmological parameters and galaxy 
bias. Starting with a fiducial cosmological model and galaxy HOD, we calculate spatial clustering 
observables on a range of length and mass scales, dynamical clustering observables that depend on 
galaxy peculiar velocities, and the galaxy-matter cross-correlation measurable by weak lensing. We 
then change one or more cosmological parameters and use y^-minimization to find the galaxy HOD 
that best reproduces the original clustering. Our parameterization of the HOD incorporates a flexible 
relation between galaxy occupation numbers and halo mass and allows spatial and velocity bias of 
galaxies within dark matter halos. Despite this flexibility, we find that changes to the HOD cannot 
mask substantial changes to the matter density Qm, the matter clustering amplitude erg, or the shape 
parameter T of the linear matter power spectrum — cosmology and bias are not degenerate. With the 
conservative assumption of 10% fractional errors, the set of observables considered here can provide 
^ 10% (Icr) constraints on as, Clm, and T, using galaxy clustering data alone. The combination 
is constrained to ~ 5%. In combination with traditional methods that focus on large-scale 
structure in the “perturbative” regime, HOD modeling can greatly amplify the cosmological power of 
galaxy redshift surveys by taking advantage of high-precision clustering measurements at small and 
intermediate scales (from sub-Mpc to ~ 20/i“^Mpc). At the same time, the inferred constraints on 
the galaxy HOD provide valuable tests of galaxy formation theory. 

Subject headings: cosmology: theory - dark matter - galaxies: formation - galaxies: halos - large-scale 
structure of universe 


1. INTRODUCTION 

From the 1970s through the early 1990s, studies of 
galaxy clustering drove much of the progress in cos¬ 
mology. Measurements of steadily improving dynamic 
range and precision demonstrated good agreement with 
the predictions of a cosmological model incorporating 
scale-invariant, Gaussian primeval fluctuations modu¬ 
lated by the transfer function expected in a universe 
dominated by cold dark matter (GDM) with flmh ~ 
0.2 (where Vim is the matter density parameter and 
h = Hq/IOO kms“^Mpc~^). The advent of multi-fiber 
galaxy redshift surveys and improved photometric in¬ 
put catalogs has dramatically improved the precision of 
clustering measurements over the last decade, beginning 
with the Las Campanas Redshift Survey (LGRS; Sheet- 
man et al. 1996) and continuing with the Two-Degree 
Field Galaxy Redshift Survey (2dFGRS; Colless et al. 
2001) and the Sloan Digital Sky Survey (SDSS; York et 
al. 2000; Abazajian et al. 2004). In parallel, numerical 
simulations and numerically tested analytic approxima¬ 
tions have turned the task of calculating non-linear dark 
matter clustering from specified initial conditions into 
an essentially solved problem. The principal obstacle to 
drawing cosmological inferences from galaxy clustering 
measurements is now the uncertainty in the relation be¬ 
tween the distribution of observable galaxies and the un¬ 
derlying distribution of dark matter, the problem known 
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as galaxy bias. Much of the cosmological progress in the 
last decade has been driven by observations that circum¬ 
vent this complication, such as cosmic microwave back¬ 
ground (CMB) anisotropies, weak gravitational lensing, 
the Lya forest, and the Type la supernova diagram. 

These new observables favor an inflationary, low-Dm, 
CDM-dominated model similar to that originally sug¬ 
gested by galaxy clustering, which in turn implies that 
the galaxies that dominate typical optically selected 
galaxy surveys must be approximately unbiased, in the 
sense that the rms galaxy count fluctuations are simi¬ 
lar to the rms dark matter density fluctuations on large 
scales (see, e.g., Lahav et al. 2002). However, observed 
galaxy clustering varies systematically with galaxy lumi¬ 
nosity, color, and spectral or morphological type (Nor- 
berg et al. 2002; Zehavi et al. 2005, and numerous ref¬ 
erences therein), and reproducing the observed galaxy 
correlation function in an inflationary CDM model re¬ 
quires that the bias of the correlation function vary with 
separation on scales below a few megaparsecs (Jenkins et 
al. 1998; Zehavi et al. 2004). The advances in observa¬ 
tional cosmology over the last few years have also raised 
the stakes for galaxy clustering studies. We are no longer 
interested in, for example, distinguishing flm ^ 0.3 from 
flm = 1; instead, we want to constrain at the few per¬ 
cent level to increase the power of tests for the nature of 
dark energy. Despite improvements in semi-analytic and 
numerical modeling of galaxy formation, it is not clear 
that these methods predict galaxy bias robustly enough 
for this kind of precision cosmology. Faced with these 
challenges, most cosmological applications of the 2dF- 
GRS and the SDSS have focused on the linear or near- 
linear regime, where generic arguments suggest that the 
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effects of galaxy bias should be relatively simple. These 
“perturbative” analyses of large-scale structure play a 
significant role in the current web of cosmological con¬ 
straints (e.g., Percival et al. 2002; Spergel et al. 2003; 
Tegmark et al. 2004b; Cole et al. 2005; Seljak et al. 2005a; 
Tegmark et al. 2006), but they are restricted to large 
scales where even these enormous surveys have limited 
statistical precision. 

In this paper we argue that recent developments in 
the theoretical description of galaxy bias allow a more 
aggressive approach to inferring cosmological constraints 
from galaxy clustering measurements. We work in the 
framework of the halo occupation distribution (HOD), 
which characterizes galaxy bias in terms of the probabil¬ 
ity distribution P{N\M) that a dark matter halo of virial 
mass M contains N galaxies of a specified type, together 
with prescriptions for the spatial and velocity bias of 
galaxies within dark matter halos (Ma & Fry 2000; Pea¬ 
cock & Smith 2000; Seljak 2000; Scoccimarro et al. 2001; 
Berlind & Weinberg 2002; Cooray & Sheth 2002). Here 
the term “halo” refers to a bound dark matter structure 
of typical overdensity pjp ~ 200, in approximate dy¬ 
namical equilibrium, which may be the individual halo 
of a single bright galaxy or the common halo of a galaxy 
group or cluster."* The flow diagram in Figure 1, adapted 
from Weinberg (2002), sketches the interplay between the 
“cosmological model” and the “physics of galaxy forma¬ 
tion” in determining observable galaxy clustering, which 
we take to include both the traditional statistics mea¬ 
sured from redshift surveys and the galaxy-matter corre¬ 
lations measured by galaxy-galaxy lensing (e.g., Fischer 
et al. 2000; Hoekstra, Yee, & Gladders 2001; Sheldon et 
al. 2004). On the left side, the cosmological model, which 
specifies the initial conditions and the energy and mat¬ 
ter contents of the universe, determines the mass func¬ 
tion, spatial correlations, and velocity correlations of the 
dark halo population. The intervening box indicates that 
the only features of the cosmological model that really 
matter in this context are and the amplitude and 
shape of the linear matter power spectrum P{k), here 
represented by erg (the rms linear matter fluctuation in 
8/i“* Mpe spheres), the inflationary spectral index n^, 
and the transfer function shape parameter F, which it¬ 
self depends on the values of flm, h, and the baryon den¬ 
sity (Bardeen et al. 1986; Hu & Sugiyama 1996). Other 
features of the cosmological model, such as the energy 
density and equation of state of the vacuum component, 
may have an important impact on other cosmological ob¬ 
servables or on the history of matter clustering, but they 
have virtually no effect on the halo population at z = 0, 
if the shape of P{k) and the present day value of erg are 
held fixed (see Zheng et al. 2002). 

On the right side of Figure 1, the box titled “galaxy 
formation physics” represents the additional processes — 
such as shock heating, radiative cooling, star formation, 
feedback, and mergers — that are essential to produc¬ 
ing galaxies and determining their masses, luminosities, 
diameters, colors, and morphologies. These physical pro¬ 
cesses operate in the background provided by the evolv- 

We have in mind the kinds of structures identified in 7V-body 
simulations by a friends-of-friends algorithm with linking length 
I ~ 0.15 — 0.2n~^P, but the precise definition of halo does not 
matter provided that one is consistent throughout all calculations. 



Fig. 1.— Interplay between the cosmological model and galaxy 
formation physics in determining observable galaxy clustering. The 
cosmological model determines the mass function and clustering of 
the dark halo population. Galaxy formation physics, operating 
within this cosmological model, determines the HOD of different 
galaxy classes. The clustering of any given class of galaxies can be 
predicted from the halo population and the HOD. In this paper 
we investigate how well one can use observations of clustering and 
galaxy-mass correlations to infer cosmological and HOD parame¬ 
ters simultaneously. 

ing halo population, so together with the cosmological 
model they determine the HODs of different classes of 
galaxies. The halo population and the HOD together de¬ 
termine galaxy clustering and galaxy-mass correlations. 
One can, of course, use hydrodynamic simulations or 
semi-analytic models to predict galaxy clustering statis¬ 
tics directly, without computing the HOD as an inter¬ 
mediate step (e.g., Kauffman et al. 1999; Benson et al. 
2000; Cen & Ostriker 2000; Pearce et al. 2001; Weinberg 
et al. 2004). However, these predictions reflect the com¬ 
bination of the cosmological model and the galaxy for¬ 
mation theory, and unless one has complete confidence 
in the latter, one cannot draw secure conclusions about 
the former. We would therefore like to know how well 
cosmological parameters can be constrained without re¬ 
lying on a detailed theory of galaxy formation, by using 
the data themselves to determine the relation between 
galaxies and dark matter. 

For this purpose, the HOD formulation of bias has 
two key strengths, both emphasized by Berlind & Wein¬ 
berg (2002). The first is the division of labor implied 
by Figure 1: galaxy formation physics influences the 
HOD, but the properties of the halo population defined 
at p/p ~ 200 are determined almost entirely by the much 
simpler physics of gravitational clustering, which can be 
modeled accurately using Y^-body simulations or numeri¬ 
cally tested analytic approximations. The second is com¬ 
pleteness: given a cosmological model and a specified 
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HOD, one can calculate any galaxy clustering statistic 
on any scale, by populating A^-body halos or by drawing 
on a steadily expanding array of analytic techniques. The 
HODs of different galaxy classes therefore encode, statis¬ 
tically, all the aspects of galaxy formation physics that 
are relevant to predictions of galaxy clustering. (We dis¬ 
cuss an important caveat to this statement below.) HOD 
modeling enables one to take full advantage of galaxy 
clustering data, with no restriction to large scales or par¬ 
ticular statistics. 

The goal of our empirical approach is to reverse the 
causal arrows in Figure 1, working backwards from the 
data to properties of the HOD and the halo population, 
and from there to conclusions about galaxy formation 
physics and fundamental aspects of cosmology. Berlind 
& Weinberg (2002) considered a fixed cosmological model 
and showed how changes to the HOD affect many of the 
traditional statistical measures of galaxy clustering, such 
as correlation functions, the group multiplicity function, 
and redshift-space distortions. They argued that the 
complementary information from different statistics and 
different spatial scales should permit an accurate empiri¬ 
cal reconstruction of the HOD for a specified cosmology, 
yielding physically informative tests of galaxy formation 
theories. Zheng et al. (2002) showed that Hm, erg, and 
the power spectrum shape have non-degenerate effects 
on the halo population, and they argued that changes 
to the galaxy HOD could not fully mask the effects of 
cosmological parameter changes. Here we complete this 
theoretical program by showing how galaxy clustering 
measurements can constrain simultaneous changes to the 
cosmological model and the galaxy HOD. Our results 
add quantitative teeth to the qualitative arguments and 
speculations in Berlind & Weinberg (2002), Zheng et al. 
(2002), and Weinberg (2002). 

We define our cosmological parameter space by the 
values of 0^, erg, and T. We concentrate mainly on 
rtm and CTg because F (or, more generally, the shape of 
the linear power spectrum) can be well constrained by 
combining the large-scale galaxy power spectrum with 
other observables like CMB anisotropy and the Lya for¬ 
est. We first compute predicted values of a number of 
galaxy clustering observables, assuming a central model 
with observationally motivated choices of cosmological 
parameters and the galaxy HOD. We then change one 
or more cosmological parameters, and we test how well 
this new cosmological model can reproduce the original 
clustering “data” given complete freedom to vary the 
HOD within a very flexible parameterization. To gain 
insight, we first explore several interesting axes within 
the (Om,crg,r) parameter space, varying parameters in¬ 
dividually or in physically motivated combinations. We 
arrive at definite numbers by assuming that each of our 
30 observables, some of which represent the same clus¬ 
tering statistic measured at multiple scales, can be mea¬ 
sured with 10% uncorrelated fractional uncertainty from 
a survey like the SDSS. This assumption seems roughly 
plausible, but the ultimate strength of cosmological con¬ 
clusions will depend on the precision and dynamic range 
of the measurements. Forecasting this precision for a 
survey like the SDSS, including the covariance of errors 
among different observables, is itself a major theoretical 
task, which we will not undertake here. For our theoreti¬ 
cal investigation, we further restrict ourselves to observ¬ 


ables for which we have reasonable analytic approxima¬ 
tions, and we suspect that our quantitative conclusions 
will in the end prove overly pessimistic because we must 
omit some observables that contain significant additional 
information. 

As discussed by Berlind & Weinberg (2002), the com¬ 
pleteness of the HOD as a description of bias rests on 
the assumption that the galaxy content of a halo of virial 
mass M is statistically independent of the halo’s larger 
scale environment. Berlind et al. (2003) show that this 
assumption accurately describes the galaxy population 
in Weinberg et al.’s (2004) hydrodynamic cosmological 
simulation. It has fundamental theoretical roots in the 
excursion set model of Bond et al. (1991), which predicts 
that the statistical features of a halo’s assembly history 
depend only on its present mass; indeed, all semi-analytic 
galaxy formation models that use statistically generated 
merger trees make this assumption implicitly. Since the 
halo mass function itself varies with environment, with 
high-mass halos absent in low-density regions, models 
that tie galaxy populations to halo masses still predict 
strong correlations between galaxy properties and the 
large-scale environment, in good agreement with obser¬ 
vations (e.g., Benson 2001; Berlind et al. 2005; Zehavi et 
al. 2005). However, while Lemson & Kauffmann (1999) 
show that Wbody halos of mass M > 10^^Mq have sim¬ 
ilar properties and formation histories in different envi¬ 
ronments, recent studies show a substantial correlation 
of halo formation redshift with large-scale overdensity 
for lower mass halos (M < IO^^’^Mq), contradicting the 
simplest form of the Bond et al. (1991) model (Gao et 
al. 2005; Harker et al. 2006; see also Sheth & Tormen 
2004). As discussed below in §5, we expect the impact of 
such dependence on the clustering statistics of mass- or 
luminosity-thresholded galaxy samples to be small, but 
given the level of precision we eventually hope to reach, 
it will probably be necessary to allow for it when fit¬ 
ting the data. Strategies for incorporating environmental 
variations into HOD modeling will require guidance from 
the next generation of hydrodynamic and semi-analytic 
models of galaxy clustering, and we reserve this task 
for future investigation. Ultimately, the assumption of 
an environment-independent HOD, or of any particular 
parameterization of environmental dependence, must be 
tested empirically by checking that the adopted model 
can consistently explain the variations of galaxy cluster¬ 
ing with large-scale environment. 

Many current efforts to constrain cosmological param¬ 
eters with galaxy clustering data make (explicitly or im¬ 
plicitly) the much stronger assumption that galaxy bias 
can be adequately described by a linear model, 6g = bSm, 
on the relevant scales. Here Sg and 6m are the galaxy 
and dark matter density contrast, respectively, and the 
linear bias factor b may depend on galaxy type but 
is assumed to be independent of scale. In this case, 
Pg{k) = b‘^Pm{k), and the shape of the galaxy power 
spectrum provides powerful cosmological constraints in 
combination with GMB data even if b is unknown (Perci- 
val et al. 2002; Tegmark et al. 2004b, 2006). Large-scale 
redshift-space distortions measure /3 = (Kaiser 

1987; Hamilton 1998; Hawkins et al. 2003), which in 
combination with clustering measurements constrains 
crgHj);® = fl^as,g/b, where as,g is the measured rms flue- 
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tuation of the galaxy density contrast at 8h~^ Mpc. The 
mass function or X-ray temperature function of galaxy 
clusters constrains a similar combination of erg and 
(White, Efstathiou, & Frenk 1993) by a route that is in¬ 
dependent of galaxy bias. At low redshift, the combina¬ 
tion of galaxy-galaxy lensing with the galaxy-galaxy cor¬ 
relation function constrains a somewhat different com¬ 
bination of these parameters, flm/b oc assuming 

a linear bias model with ^gg(r) = b^g^{r) = 

The triangle shape dependence of the reduced galaxy 
bispectrum can constrain b directly (Fry 1994; Verde et 
al. 2002), in this case assuming a quadratic bias model 
Sg = b6m + + const, to relate the galaxy and mass 

density fields at second order (Fry 1994; Juszkiewicz et 
al. 1995). 

General theoretical arguments suggest that linear bias 
should be a good approximation for the power spectrum 
on sufficiently large scales, provided that the efficiency 
of galaxy formation is determined by the local (r < few 
Mpc) environment (Coles 1993; Fry & Gaztahaga 1993; 
Weinberg 1995; Mann, Peacock, & Heavens 1998; Scher- 
rer & Weinberg 1998; Narayanan, Berlind, & Weinberg 
2000). Numerical experiments using such “local” bias 
models also show that the “b” multiplying the power 
spectrum amplitude should be similar to the “b” affecting 
redshift-space distortions (Berlind, Narayanan, & Wein¬ 
berg 2001). Thus, the combination of these perturbative 
galaxy clustering analyses can in principle yield separate 
constraints on T (measured directly from the power spec¬ 
trum shape), flm, and ag, with the degeneracy of flm 
and (Tg broken either by the bispectrum analysis or by 
the different flm dependence of the galaxy-galaxy lens¬ 
ing constraint and the /3 or cluster normalization con¬ 
straints. The applications of these techniques to 2dF- 
GRS and SDSS data show an impressive degree of inter¬ 
nal consistency and good agreement with external cos¬ 
mological constraints (Percival et al. 2002; Verde et al. 
2002; Spergel et al. 2003; Tegmark et al. 2004b; Cole et 
al. 2005; Seljak et al. 2005a; Tegmark et al. 2006). How¬ 
ever, reliance on the linear or quadratic bias approxi¬ 
mations restricts these analyses to large scales, and it is 
not clear that these approximations hold consistently at 
the level of accuracy desired for further improvements 
(better than 10%, say), since scatter in the relation be¬ 
tween galaxy and mass densities can have different effects 
on different statistics and at different scales (Pen 1998; 
Dekel & Lahav 1999). The HOD approach to galaxy clus¬ 
tering analysis substitutes a much more general model of 
galaxy bias, and it makes use of high-precision measure¬ 
ments from small and intermediate scales (from sub-Mpc 
to ~ 20h~^ Mpc) in addition to the lower precision mea¬ 
surements in the perturbative regime. Furthermore, re¬ 
gardless of the strength of the cosmological constraints, 
the HOD parameters derived from the data themselves 
provide detailed tests of galaxy formation models. 

The program of HOD-based interpretation of observed 
galaxy clustering is, in fact, well underway. Jing & 
Borner (1998) and Jing, Borner, & Suto (2002) used 
HOD-type bias models to interpret the correlation func¬ 
tions and pairwise velocity dispersions measured from 
the LGRS and the Point-Source Catalog Redshift Sur¬ 
vey (PSCz; Saunders et al. 2000). Peacock & Smith 
(2000), Marinoni & Hudson (2002), and Kochanek et al. 
(2003) used the group multiplicity function to constrain 


the galaxy occupations of high-mass halos. Guzik & Sel¬ 
jak (2002), Seljak et al. (2005a), and Mandelbaum et al. 
(2006) applied HOD modeling to galaxy-galaxy lensing 
measurements from the SDSS. Zehavi et al. (2004), an¬ 
alyzing a volume-limited sample of bright {Mr < —21) 
SDSS galaxies, showed that HOD models naturally ex¬ 
plain the observed deviation from a power-law correla¬ 
tion function. Zehavi et al. (2005) measured the lumi¬ 
nosity and color dependence of the SDSS galaxy corre¬ 
lation function and used it to infer the luminosity and 
color dependence of the galaxy HOD, finding results in 
good qualitative agreement with theoretical predictions 
(Zheng et al. 2005). Magliocchetti & Porciani (2003) 
used a similar approach to investigate the halo occupa¬ 
tions of early- and late-type galaxies in the 2dFGRS, 
while Collister & Lahav (2005) derived HODs of red and 
blue 2dFGRS galaxies from the group catalog of Fke 
et al. (2004), testing the consistency of their result by 
comparing predicted and observed correlation functions. 
Abazajian et al. (2005) used the SDSS measurements for 
Mr < —21 galaxies in conjunction with CMB anisotropy 
data to infer simultaneous constraints on HOD and cos¬ 
mological parameters. Tinker et al. (2005) used HOD 
modeling of Zehavi et al.’s (2005) clustering measure¬ 
ments to predict cluster mass-to-light {M/L) ratios, and 
they inferred constraints on by comparing to pub¬ 

lished M/L measurements. HOD models have been ap¬ 
plied to the interpretation of high-redshift clustering by 
Bullock et al. (2002), Moustakas & Somerville (2002), 
Yan, Madgwick, & White (2003), Zheng (2004), Ouchi 
et al. (2005), Lee et al. (2006), Coil et al. (2006), and 
Cooray (2006). Finally, van den Bosch, Mo, & Yang 
(2003a) have initiated a comprehensive program simi¬ 
lar to the one described here, based on the closely re¬ 
lated conditional luminosity function (CLF) formalism 
(see also van den Bosch et al. 2003b, 2004, 2005, 2006; 
Yang, Mo, & van den Bosch 2003; Yang et al. 2004, 2005; 
Mo et al. 2004; Wang et al. 2004). We discuss the sim¬ 
ilarities and differences between the HOD and CLF ap¬ 
proaches in § 6. 

In the next section we define our class of cosmological 
models, list the analytic approximations we use for prop¬ 
erties of the halo population, and describe our flexible pa¬ 
rameterization of the galaxy HOD. In § 3 we describe the 
methods and approximations that we adopt for analytic 
calculation of galaxy clustering observables. Section 4 
presents our main results, showing the ability of com¬ 
plementary galaxy clustering measurements to constrain 
the galaxy HOD and cosmological parameters simulta¬ 
neously. Section 5 discusses the issue of environmental 
variations of the HOD. Section 6 summarizes our Hnd- 
ings and discusses the overall prospects for application 
of this approach. A reader familiar with the field who 
wants an overview of our main results can read the Hrst 
paragraph of §2.1, skim § 2.2 to understand our HOD 
parameterization, then skip to § 4, paying particular at¬ 
tention to § 4.2, § 4.5, and § 4.6 and Figures 5, 14, and 
15. 

2. COSMOLOGICAL MODEL AND HOD 
PARAMETERIZATION 

2.1. Cosmological Parameters and Halo Properties 

Throughout this paper we adopt spatially flat ACDM 
cosmological models with Gaussian initial density flue- 
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tuations. The cosmological model is defined by the mass 
density parameter Qm and the mass fluctuation power 
spectrum P{k) cx T), where k is the wavenum¬ 

ber, Us is the spectral index of the inflationary power 
spectrum, and T(fc; T) is the transfer function with shape 
parameter T. We parameterize each cosmological model 
in this paper by flmi erg, Ug, and T, where erg is the 
rms fluctuation of the linear density field filtered with a 
top-hat filter of radius 8/i“^Mpc. We adopt the param¬ 
eterization of Efstathiou, Bond, & White (1992) for the 
transfer function T(k; T), which approximates the evolu¬ 
tion of adiabatic primordial fluctuations. This formula¬ 
tion suffices for our purpose of investigating how galaxy 
clustering data constrain HOD and cosmological param¬ 
eters, with a single parameter T encoding the combined 
effects of the matter-radiation transition (at the scale 
A oc ^Imh) and the suppression of fluctuation growth in 
baryon or massive neutrino components. Observational 
analyses should use a more accurate transfer function, 
e.g., from CMBFAST (Seljak & Zaldarriaga 1996), but 
this would not change the sensitivity to the shape and 
amplitude of P(k). It is important to note that we treat 
Dm and T as independently variable quantities; we do not 
automatically change the power spectrum shape when we 
change 0^. In the context of ACDM, our approach im¬ 
plicitly assumes that a change in is compensated by 
a change in h or some other parameter to keep the power 
spectrum at its empirically constrained shape. 

Our calculations of galaxy clustering statistics rely on 
analytic descriptions of the density profiles, mass func¬ 
tion, spatial clustering, and velocity statistics of dark 
matter halos. We draw on the extensive literature that 
presents numerically tested analytic approximations or 
numerically calibrated fitting formulae for these quanti¬ 
ties. 

The density profile of a dark matter halo of mass M is 
assumed to have the Navarro-Frenk-White (NFW) form 
(Navarro, Frenk, & White 1995, 1996, 1997), 


Pm{r,M) = 


(r/rs)(l + r/rs)'^’ 


where the characteristic radius is related to the virial 
radius i?vir of the halo through the concentration param¬ 
eter c = i?vir/?"s, and 


Ps — Po- 


3 ln(l-I-c) — c/(l-I-c) 


We define the virial radius by the condition Jvir = 200, 
independent of cosmology, where dvir is the average mass 
density of the halo within i?vir in units of the mean mat¬ 
ter density po- The concentration parameter c is a func¬ 
tion of halo mass. Here we adopt the relation given by 
Bullock et al. (2001) after modifying it according to our 
constant-5vir halo definition. 


7 M 

c(M) = « ^ 


where Cg = 11, /3c = —0.13, and M* is the nonlinear 


The mass function of halos can be expressed as 
dM M ^ ' dM 


where n{M) is the space density of halos with mass 
higher than M, a{M) is the rms fluctuation of the 
mass overdensity at a mass scale M, v = 5c/cr{M) with 
5c ~ 1.686 is the threshold density contrast for collapse, 
and f{i') is a dimensionless function. The mass M* in 
equation (3) is defined by the condition = Sc (i.e., 

V = 1). The function /(jz) can be approximately derived 
using the Press-Schechter (1974) or excursion set (Bond 
et al. 1991) formalism, or it can be calibrated by fitting 
the results of A^-body simulations (e.g., Sheth & Tor- 
men 1999; Jenkins et al. 2001). Here we use the fitting 
formula of Sheth & Tormen (1999), 

f{v)=A^l^ l + :/exp^-^^, (5) 

where A = 0.3222, a = 0.707, and p = 0.3. This formula 
is consistent with the ellipsoidal collapse model (Sheth, 
Mo, & Tormen 2001), and it provides an excellent match 
to results from simulations (see, e.g., Reed et al. 2003). 

At large scales, the halo two-point correlation function 
^hh(c) is biased relative to the matter two-point correla¬ 
tion function ^mm(?") by a mass-dependent factor bf(M). 
We adopt Sheth et al.’s 2001 formula for bh{M), based 
on the ellipsoidal collapse model, 

bh{M) = l-\ —ri— Va(ai^^)-I-^/a6(a^^)^“° 

^ya^c _ 

__ 

{av'^Y + ^(1 ~ c)(l — c/2) 

where a = 0.707, b = 0.5, c = 0.6, and v and 5c have 
the same meaning as in equation (4). Our calculation 
of the galaxy 2-point correlation function, discussed in 
§ 3.1, also incorporates scale dependence of halo bias. 
Recent numerical work (Seljak & Warren 2004; Tinker 
et al. 2005) suggests that the bias factors given by equa¬ 
tion (6) are systematically too high. Since we use the 
same approximations to calculate the “observations” of 
our central model and the “predictions” of others, we do 
not expect moderate errors in bh{M) to alter our conclu¬ 
sions. However, comparisons to real observational data 
should use the most accurate available bias factors. 

The rms three-dimensional (3D) space velocity of halos 
of mass M can be approximated as 

ah{M) = (7) 

where 

1 

a‘^(M) = — dk k‘^+‘^^P{k)W‘^[kR{M)], (8) 

27r Jg 

W[kR{M)] is the Fourier transform of the top-hat filter 
on mass scale M, and the smoothing radius R is deter¬ 
mined through M = 4TrpgR^/3 (Sheth & Diaferio 2001). 
The one-dimensional (ID) rms velocity is simply the 3D 
one divided by -s/S- The rms ID velocity difference of 
halos of masses Mi and M 2 at a separation r is 

A72|r) = a'l iY){Mi)+al^ iD{M2)+'^{Mi,M2\r), 

(9) 

where 4'(Mi,M2|r) represents the correlation term. In 
our calculations (see §3.2), we ignore the correlation 
term, which is an acceptable simplification (Sheth & Di¬ 
aferio 2001). We assume that the ID velocity disper¬ 
sion of dark matter within a halo of mass M has the 
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value (G'M/2i?vir)^^^ expected for a singular isothermal 
sphere, independent of radius. Moderate departures from 
this assumption can be accommodated by the velocity 
bias parameter ay, described below. 

2.2. HOD Parameterization 

The HOD of a specified class of galaxies is defined by 
the probability distribution P{N\M) that a halo of mass 
M hosts N galaxies, together with prescriptions for the 
spatial and velocity distributions of galaxies within ha¬ 
los relative to those of the dark matter particles. For 
our purpose in this paper, we want a flexible parame¬ 
terization of the HOD that could describe the relevant 
features of any physically reasonable model. We want the 
constraints on the HOD to be imposed by the clustering 
measurements rather than by our parameterization, and 
we want to give each cosmological model its best shot at 
reproducing these measurements before rejecting it. 

Conceptually, it is useful to separate P{N\M) 
into a mean occupation function {N{M)), defined 
by {N{M)) = J2n ^ ^ distribution 

P{N\{N)) at fixed halo mass. There have been numerous 
investigations of the predictions of semi-analytic mod¬ 
els, hydrodynamic simulations, and high-resolution N- 
body simulations for P{N\M) (e.g., Kauffmann, Nusser, 
& Steinmetz 1997; Benson et al. 2000; Seljak 2000; Scoc- 
cimarro et al. 2001; Sheth & Diaferio 2001; White, Hern- 
quist, & Springel 2001; Yoshikawa et al. 2001; Cooray & 
Sheth 2002; Guzik & Seljak 2002; Scranton 2003; Berlind 
et al. 2003; Kravtsov et al. 2004; Zheng et al. 2005; the 
last three are the most extensive and the most relevant 
to our considerations here). The predicted HOD depends 
on the defining characteristics of the galaxy sample, such 
as luminosity, color, and morphological type. For the 
most sensitive cosmological investigations, it seems best 
to work with volume-limited samples defined by luminos¬ 
ity thresholds, which are the largest homogeneous sam¬ 
ples (i.e., with the same distribution of galaxy types at 
all distances) that one can construct from an apparent 
magnitude-limited survey. In this case, {N{M)) should 
be a more or less monotonic function of halo mass, since 
more massive halos contain, on average, either more 
galaxies or brighter galaxies. Below some cutoff mass, 
{N{M)) should drop to zero, as low-mass halos do not 
contain enough cooled baryons to make a galaxy above 
the luminosity threshold. Instead of a luminosity thresh¬ 
old, one can impose a threshold in stellar mass, esti¬ 
mated for each galaxy from its luminosity and its color 
or spectral energy distribution (e.g.. Bell & de Jong 2001; 
Kauffman et al. 2003). This approach should reduce the 
scatter between the galaxy observable and the host halo 
mass by reducing the impact of stellar population vari¬ 
ations, thus sharpening the cutoff in {N{M)) and min¬ 
imizing its environmental dependence. However, if one 
starts from an apparent magnitude-limited survey, then 
the mass-thresholded samples one can create are smaller 
than the luminosity-thresholded samples because of the 
range of galaxy colors. 

While we do not want to impose strong theoretical pri¬ 
ors on our HOD parameterization, we have decided to in¬ 
corporate one of the most generic and useful results that 
comes from the galaxy formation papers cited above: the 
distinction between central and satellite galaxies. In hy¬ 
drodynamic simulations and semi-analytic models, most 


halos contain a galaxy near the center of mass, moving 
close to the center-of-mass velocity, which is usually more 
massive and older than any other galaxies in the halo 
(Berlind et al. 2003). Kravtsov et al. (2004) show that 
the central-satellite distinction naturally explains one of 
the quantitatively important properties of the HOD, the 
sub-Poisson width of P{N\{N)) at small mean occupa¬ 
tion numbers. For satellite dark matter subhalos, the 
fluctuations about the mean occupation are close to Pois¬ 
son, but the number of central galaxies is by definition ei¬ 
ther zero or one, thus obeying nearest-integer (Bernoulli) 
statistics. Fluctuations in the full P{N\{N)) are substan¬ 
tially sub-Poisson as long as the central galaxy makes a 
significant contribution to {N{M)). Zheng et al. (2005) 
show that this argument carries over to galaxies in hy¬ 
drodynamic simulations and semi-analytic models. 

Motivated by these results, we model the mean occu¬ 
pation function for galaxies above a luminosity threshold 
as 

{N{M)) = i [l + erf [1 + 5(M)], 

( 10 ) 

where the last set of square brackets represents the sum 
of one central galaxy and S{M) satellites and the hrst 
provides a smooth cutoff near M^nin- The form of the 
cutoff prohle corresponds to a log-normal distribution of 
central galaxy luminosity at hxed halo mass (see Zheng 
et al. 2005). We apply the same cutoff profile to satellite 
galaxy numbers, although in practice these are already 
well below unity close to Aimin- To model S{M) in a flex¬ 
ible way, we take the free parameters to be the values of 
log S'(M) at m fixed values of logM and dehne the con¬ 
tinuous function by a cubic spline that passes through 
these TO values. We set to = 5 and choose values of 
logS'(M) at \og{M/h~^M q) = 11, 12, 13, 14, and 15 for 
the spline ht. Halos more massive than Mq are 

too rare to have much impact on the clustering statis¬ 
tics we consider here. Theoretical predictions of S{M) 
can be accurately described by a truncated power-law 
parameterized by a slope, amplitude, and low-mass cut¬ 
off (Kravtsov et al. 2004; Zheng et al. 2005), but our 
five-parameter model allows for curvature or inflections 
that depart from these theoretical predictions. In addi¬ 
tion to the S{M) values, the adjustable parameters of 
equation (10) are CTiogM and M^-m- We always assume 
that the uncertainty in the galaxy mean density is negli¬ 
gible in comparison to the uncertainty in the clustering 
measurements. Therefore, once we have specified S{M) 
and O’log M, we choose the value of that yields the 
correct mean number density, rather than treating it as 
a free parameter. (We could pick any of the parameters 
for this treatment, but the mean density is most sensitive 
to Mmin-) 

For the clustering measures and analytic approxima¬ 
tions used in this paper, the only property of P{N\{N)) 
that we need to specify is the second factorial moment 
{N{M)[N(M) — 1]), which determines the mean num¬ 
ber of pairs per halo. The number of central galaxies is 
zero or one with a relative probability determined by the 
mean {Ncen{M)). For satellite galaxies, we introduce a 
width parameter lo defined by 

(fVsat(M)[Afsat(M) - 1]) = w(fV,at(M))2. (11) 

The value uj = 1 corresponds to a Poisson distribution 
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of satellite numbers NgatiM) at fixed halo mass, con¬ 
sistent with theoretical predictions. The parameter u 
allows the satellite probability to be broader or narrower 
than Poisson. The full probability distribution becomes 
narrow at low {N{M)) in any case because of central 
galaxies. The value of u determines how quickly the 
distribution approaches a Poisson-like width, thus influ¬ 
encing the number of galaxy pairs in halos with mean 
occupation numbers ~ 1—several. 

In the hydrodynamic simulation examined by Berlind 
et al. (2003), satellite galaxies have a radial profile sim¬ 
ilar to that of the dark matter. However, the physical 
processes that determine the satellite profile are complex 
(see, e.g., Nagai & Kravtsov 2005; Zentner et al. 2005), 
so there is no reason to expect the profiles to be identi¬ 
cal. Here we assume that satellite galaxies in a halo still 
follow an NFW profile, but we allow the galaxy concen¬ 
tration parameter Cg to differ from the underlying dark 
matter concentration parameter c. We characterize the 
relative spatial distributions of satellite galaxies and dark 
matter by A logco = logCg—logco, where a positive (neg¬ 
ative) A log Co means that galaxies are more (less) cen¬ 
trally concentrated than dark matter. We assume that 
central galaxies always reside at the halo center of mass. 

We also allow the velocity dispersion of satellites to 
differ from the velocity dispersion of dark matter by an 
average factor a^, which we assume to be independent 
of mass. Central galaxies are assumed to move at the 
halo center-of-mass velocity. Simulations suggest that 
av ~ but the results are not entirely consistent on this 
point (e.g., Berlind et al. 2003; Kravtsov et al. 2004), 
and in any event we want this dynamical question to be 
settled by observations rather than by theory. Tinker et 
al. (2006) discuss the influence of on redshift-space 
distortions in detail. 

Altogether we have 10 free parameters in our HOD 
prescription — seven in {N{M)) (Mmin, criogM, five 
spline points), one (w) for the second moment of 
P{N\{Nsa.tiM))), one (Alogco) for the spatial bias 
within halos, and one («„) for the velocity bias within 
halos. One of these (which in practice we take to be 
Mynin) is determined by matching the mean galaxy space 
density once the values of other parameters are specified. 
The remaining nine are varied to yield the best match 
to the observables described in §3. The parameters w, 
A log Co, and could in principle depend on halo mass, 
but the variation would have to be quite strong over a 
fairly narrow mass range to have an effect that could not 
be adequately described by an appropriate average value. 

The main theoretical preconception that is built into 
our parameterization is the existence of central galax¬ 
ies. The central location of these galaxies makes a small 
but not negligible difference to the two-point correlation 
function (see Berlind & Weinberg 2002), although this 
could probably be compensated by Alogco. However, 
since {Nsat{M)) in practice turns out to be small near 
the cutoff mass Mmin, the central galaxy parameteriza¬ 
tion necessarily leads to sub-Poisson P{N\{N)) fluctua¬ 
tions at low (A(M)). Since this aspect of galaxy forma¬ 
tion physics is well grounded, we think it is reasonable 
to adopt it, and it allows us to represent a wider range of 
physically realistic HODs with our finite parameter set. 
However, even if we do not make the central-satellite 
distinction and instead use a functional form like that of 


Scranton (2003) to describe the width of P{N\{N)) as a 
function of halo mass, we reach similar conclusions about 
the ability to break degeneracies between cosmology and 
galaxy bias. In essence, the observations drive us to sub- 
Poisson fluctuations at low halo masses whether or not 
we impose this trend as a theoretical expectation (see 
Benson et al. 2000; Peacock & Smith 2000; Berlind & 
Weinberg 2002). 

3. ANALYTIC CALCULATION OF OBSERVABLE 
QUANTITIES 

The observable quantities considered in this paper fall 
into three general categories. First are quantities re¬ 
lated to the real-space clustering of galaxies: the rms 
galaxy number density contrast, the galaxy two-point 
correlation function, the reduced galaxy bispectrum, and 
the group multiplicity function. Second are dynami¬ 
cal quantities that depend on galaxy peculiar veloci¬ 
ties: the large-scale redshift-space distortion parameter 
P, the galaxy pairwise velocity dispersion, and the es¬ 
timated virial mass of galaxy groups as a function of 
group richness. Finally, we include the galaxy-matter 
two-point cross-correlation function, which can be mea¬ 
sured by galaxy-galaxy leasing. Most of these quantities 
are continuous functions of scale or group multiplicity. 
To obtain a discrete set of observables, we sample these 
functions at a discrete set of separations or multiplicities, 
listed in §4. We space our samplings widely so that we 
can treat the statistical errors in each measurement as 
approximately uncorrelated. 

In addition to these observables, we impose the require¬ 
ment that the HOD model reproduce the correct value 
of the sample’s mean space density, implying 

pOO 7 

n,=l dM^{N{M)). (12) 

In effect, we treat fig as an observable with negligible 
statistical error. 

In this paper we limit ourselves to observables for 
which we have some reasonable analytic approximation. 
In many cases, these approximations are not accurate 
enough for application to high-precision clustering mea¬ 
surements from large surveys like the SDSS or 2dFGRS 
— or at least their accuracy has not been extensively 
tested with numerical simulations. However, we use the 
same approximations to calculate “observed” quantities 
for our central cosmological model and HOD and “pre¬ 
dicted” quantities for altered models. Provided that 
our approximations realistically capture the variations 
of the observables with variations in the cosmological 
and HOD parameters, they are adequate for our present 
purposes of forecasting the ability of clustering measure¬ 
ments to break the degeneracy between cosmology and 
galaxy bias. Observational applications of these tech¬ 
niques will require further testing and development of 
analytic approximations. They could also benefit from 
analytic approximations for other statistics, or from nu¬ 
merical studies that provide accurate interpolation for¬ 
mulae in the relevant ranges of cosmological and HOD 
parameters. 

3.1. Spatial Clustering 

At large scales, the galaxy two-point correlation func¬ 
tion Cgg(i’) is biased relative to the matter correlation 
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function by (the square of) a number-weighted average 
of the halo bias factor, 

1 dn 

= dM-(iV(M))MM) . (13) 

We use equations (4) and (5) for dn/dM and equa¬ 
tion (6) for bh{M). The fractional rms fluctuation (Tg{R) 
of galaxy number counts in spheres of radius R can be 
expressed as an integral over $gg(r), so at large R it is 
related to the rms mass fluctuation a{R) by 

o'g(^) = bga{R) . (14) 

At R = 8h~^ Mpc, the smaller of the two scales we use 
for <7g{R), the approximation (14) is good but not perfect 
(see Tinker et al. 2006). 

More generally, Cgg(''’) can be decomposed into one- 
halo and two-halo terms, Cgg(c) and ^gg(r), which re¬ 
spectively represent pairs of galaxies residing in the same 
halo and in separate halos. The one-halo term dominates 
at small separations and the two-halo term at large sep¬ 
arations, with the transition occurring near the virial 
diameter of large halos. Noting that the total num¬ 
ber of galaxy pairs [cx 1 -I- Cgg(^)] is simply the sum 
of the number of pairs from single halos [oc 1 -b Cgg(c)] 
and that from different halos [oc 1 -b C|g(^)]) we have 

4gW = [i + 4"W] + 4'W. 

The one-halo term Cgg('c) can be exactly computed in 
real space through (Berlind & Weinberg 2002) 


1 + CM = ; 


Jo 

1 


2i?vir(M) 


dM 


F' 


dn {N{M)[N{M)-l]) 


dM 


2i?v 


(15) 


where {N{M)[N{M) — l])/2 is the average number of 
pairs in a halo of mass M and F(r/2i?vir) is the cumu¬ 
lative radial distribution of galaxy pairs, i.e., the aver¬ 
age fraction of galaxy pairs in a halo of mass M (virial 
radius i?vir) that have separation less than r. The dif¬ 
ferential function F'{x) is determined by the proHle of 
the galaxy distribution within the halo. We assume that 
there is always a galaxy located at the center of the 
halo, and others are regarded as satellite galaxies. With 
this assumption, F'{x) is the pair-number weighted av¬ 
erage of the central-satellite pair distribution F^^{x) and 
the satellite-satellite pair distribution F!.^{x) (see, e.g., 
Berlind & Weinberg 2002; Yang, Mo, & van den Bosch 
2003), 


The central-satellite galaxy pair distribution F^g{x) is 
just the normalized radial distribution of galaxies within 
halos, and the satellite-satellite galaxy pair distribution 
Fgg(a;) can be derived through the convolution of the 
galaxy distribution profile with itself (Sheth & Diafe- 
rio 2001). In the case of an NFW proHle, as used in 
this paper, {x) can be analytically expressed (see Ap¬ 
pendix A; also see Sheth & Diaferio 2001). 


At large scales, the two-halo term becomes Cgg(^) = 
5g^mm(?’), where bg is the galaxy number-weighted halo 
bias factor of equation (13). At smaller scales, one must 
account for the convolution of galaxy proHles over the 
Hnite size of halos, scale dependence of the halo bias fac¬ 
tor, and the fact that separate halos by definition do 
not interpenetrate. Because convolutions become multi¬ 
plications in Fourier space, it is convenient to calculate 
the two-halo component of the galaxy power spectrum 
Pgg{k) and Fourier transform it to get Cgg(^) (Seljak 
2000; Scoccimarro et al. 2001). Our full method of cal¬ 
culating the two-halo term is described in Appendix B 
of Tinker et al. (2005). Our adopted treatment of halo 
exclusion is referred to there as the “spherical” exclusion 
method, and we use the Sheth, Mo, & Tormen (2001) 
formula for bh{M) instead of the modified formula pro¬ 
posed by Tinker et al. (2005). We have used this method 
in earlier papers (Zehavi et al. 2004, 2005; Zheng 2004). 

Galaxies within the same dark matter halo constitute 
a galaxy group. The cumulative multiplicity function of 
galaxy groups, which is the number density of galaxy 
groups that have N or more members, is 

dn 

ngrp(> N) = J2 dM—Pm)- (17) 

In this work we make the simplifying approximation 
that P{i\M) is a nearest-integer distribution, where N 
is one of the two integers bracketing {N{M)), and the 
relative probability is determined by having the right 
mean. The space density of galaxy groups with multi¬ 
plicity > {N{M)) is approximately equal to the space 
density of halos with mass > M, so given a cosmologi¬ 
cal model that specifies n{M) one can effectively “read 
off” the high end of the mean occupation function (Pea¬ 
cock & Smith 2000; Berlind & Weinberg 2002; Marinoni 
& Hudson 2002; Kochanek et al. 2003). A full observa¬ 
tional implementation should take into account both the 
width of P{N\{N)) and the scatter between estimated 
and true multiplicities introduced by the group finding 
algorithm. We consider galaxy groups with five or more 
members, so the approximation (17) is fairly good. 

3.2. Dynamical Measurements 

In redshift space, the peculiar velocities of galaxies dis¬ 
tort the galaxy power spectrum and correlation function, 
producing anisotropy in which the line of sight is a pre¬ 
ferred direction. By measuring this distortion on large 
scales, where it is caused by coherent flows into over- 
dense regions and out from underdense regions, one can 
infer a combination /3 = /bg of the density parame¬ 
ter and the large-scale galaxy bias factor (Kaiser 1987). 
We adopt P as the first of our dynamically sensitive ob¬ 
servables. The task of inferring P from clustering mea¬ 
surements in the context of HOD models is discussed in 
detail by Tinker et al. (2006). 

The next dynamical measurement we consider is 
the pairwise radial velocity dispersion of galaxies, 
defined as cr^(r) = {V 12 ) — ( 1 ^ 12 )^, where V 12 = 

— (vi — V 2 ) • (ri — r 2 )/|ri — r 2 | is the radial (inward) 
velocity of a galaxy pair, and (i = 1, 2) are veloci¬ 
ties and positions of the two galaxies, and the average is 
over all galaxy pairs with separations near r = |ri — r 2 |. 
The pairwise dispersion can be inferred by modeling the 
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redshift-space correlation function (e.g., Bean et al. 1983; 
Davis & Peebles 1983). 

Similar to the galaxy correlation function, on small 
scales, cr^(r’) is dominated by a one-halo term, which can 
be computed through (Berlind & Weinberg 2002) 


cr 


2 

D,lh 


(r) 


1 {N{M)[N{M)-1]) 

27rr2n24g(r) Jo dM 2 


2i?v 


-F' 



(18) 


where should be understood as the pair-number 

weighted average dispersion of the relative radial veloc¬ 
ities between pairs of galaxies in the same halo. If we 
assume that the velocity distribution of satellite galaxies 
is isotropic and isothermal and that the central galaxy 
is at rest with respect to the center-of-mass of the par¬ 
ent halo, then the pair-weighted average dispersion can 
be separated into central-satellite and satellite-satellite 
terms as in equation (16), 


{N{M)[N{M) - 1]) 


F'alJM) = {NUM))F'^,alal{M) 


(iVsat(M)[iV,at(M) - 1]) 


FiJalaliM) . 


(19) 


Here = GM/2i?vir is the ID velocity dispersion 

of dark matter within the halo, ay is the velocity bias 
factor, and other quantities have the same meaning as in 
equation (16). Note that there is a factor of 2 for satellite- 
satellite velocity dispersions, since both members of the 
pair are moving with respect to the halo center-of-mass. 

On large scales, the two galaxies of each pair come 
from different halos, so the relative motions of halos also 
contribute to the galaxy pairwise dispersion. We ignore 
the halo velocity correlation (d> in eq. [9]), which is an 
acceptable simplification (Sheth & Diaferio 2001). Fol¬ 
lowing the argument of Berlind & Weinberg (2002) but 
taking into account the separation of central and satellite 
galaxies, the galaxy pairwise radial velocity dispersion on 
large scales can be expressed as 




^g[l+4gW]io dMi 

dn 

X J dM2-^jj^[l + ^hh{r]Mi,M2)] 

X [(fV(Mi)) {NiM2))al,ooir; Mi, M 2 ) 

+ (7V,at(Mi))(iV(M2))a2a^(Mi) 

+ (iV(Mi))(7V,at(M2))a^a^(M2)], (20) 


where ^hhir; Mi, M 2 ) is the two-point correlation func¬ 
tion between halos of mass Mi and M 2 at separation r. 
The three terms inside the last set of square brackets rep¬ 
resent contributions by the pairwise velocity dispersion of 
Ml and M 2 halos, by the dispersion of satellite galaxies in 
Ml halos, and by the dispersion of satellite galaxies in M 2 
halos. The halo-halo velocity dispersion is imprinted on 
all {N{Mi)){N{M 2 )) galaxy pairs. However, the veloc¬ 
ity dispersion of satellite galaxies in, say. Mi halos only 
contributes to the dispersion of the {Nsat{Mi)){N{M 2 )) 
pairs involving those satellites. The halo-halo velocity 
dispersion can be computed through equation (9), and on 
large scales we have &h(r; Mi, M 2 ) = &?i(Mi)&?i(M 2 )Cmm 


and Cgg(^) = With these substitutions, equa¬ 

tion (20) reduces to 
2 

0'«,2h(^) = 1 , y + ?gg(’")('^MD)jV6] 

^ SggHl 

20-2 

+ T-j^[{<Tl)N.+U,{r){al)M.bWl) 

^ sggHl 

where the number- and bias-weighted averages {)n and 
0 Nb are defined as 

{f)N = ^ r dM-^{N{M))f{M), (22) 

Ug Jo dM 
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{f)Nb=^ dM^{N{M))bh{M)f{M), (23) 

Ugbg Jo dM 

and ()jv^ and {)Nsb are defined by replacing {N{M)) with 
(A^sat(M)) in equations (22) and (23), respectively. Since 
we do not have a good approximation for (r) at inter¬ 
mediate scales, we only consider pairwise dispersion mea¬ 
surements at small scales, where equation (18) applies, 
or on large scales, where equation (21) applies. Ulti¬ 
mately, the constraints obtainable from anisotropies in 
redshift space on small, intermediate, and large scales 
can be inferred using the methods of Tinker et al. (2006) 
and Tinker (2007). 

We also consider dynamical measurements of the av¬ 
erage virial masses of galaxy groups of fixed multiplicity, 
(Mvir(lV)). These can be inferred from galaxy positions 
and velocities using the kind of estimators described by 
Heisler, Tremaine, & Bahcall (1985). These estimators 
are affected by velocity bias, with the estimated mass 
scaling as a^. For a given halo mass function and HOD, 
we therefore calculate this quantity as 


(Mvir(lV)) 


2 /o“ dM{dn/dM)P{N\M)M 
dM{dn/dM)P{N\M) 


(24) 


As in our multiplicity function calculation, we approxi¬ 
mate P{N\{N)) by a nearest-integer distribution. With 
equation (24), we implicitly assume that the mass esti¬ 
mators are applied only to satellite galaxies, or else that 
the reduced dispersion due to the central galaxy is prop¬ 
erly taken into account. With large enough galaxy num¬ 
bers, more sophisticated mass estimators can circumvent 
the effects of velocity bias (e.g., Carlberg, Yee, & Elling- 
son 1997). In addition, average group masses can be 
estimated using gravitational lensing or X-ray properties 
instead of galaxy velocities. We do not consider these 
kinds of ai,-free mass estimates in our analysis, but it 
is clear that they would add power to our overall set of 
constraints if they could be implemented robustly. 


3.3. Galaxy-Mass Clustering 

The galaxy-mass cross-correlation can be probed by 
galaxy-galaxy lensing (e.g., Fischer et al. 2000; Hoek- 
stra, Yee, & Gladders 2001; Sheldon et al. 2004; Seljak 
et al. 2005a). In practice, what is usually inferred from 
galaxy-galaxy lensing at low redshift is the product of the 
mean matter density flm and the galaxy-mass two-point 
cross-correlation function ■Cgm(?"), so we adopt values of 
Dm^gm(r) as the observables in our analysis. 

The galaxy-mass two-point cross-correlation function 
can also be decomposed into one-halo and two-halo 
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terms. For simplicity, we only consider small-scale cross¬ 
correlation, where the one-halo term dominates. Similar 
to the one-halo term of the galaxy two-point correlation 
function, the one-halo term ^gm(r) of the galaxy-mass 
two-point cross-correlation function is computed in real 
space by counting galaxy and mass particle pairs, 

1 -|- f (r) = - 

gm Airr'^ngQrnPc 

X [ dM—{N(M))M - \ — -F' ( V,25) 

where pc is the critical density. The function Fg^{x) 
is the cumulative radial distribution of galaxy and mass 
particle pairs (not galaxy-galaxy pairs as in eq. [15]), 
which is determined by the distribution profiles of galax¬ 
ies and dark matter within the halo. When the effect of 
the central galaxy is taken into account, we have (similar 
to eq. [16]) 

(iV(M))F'^(a:) = FU^) + (iV,at(M))F;„(x), (26) 

where F^^ (Fsm) is the cumulative distribution of pairs 
made of the central galaxy (satellite galaxies) and mass 
particles. The function F^^{x) is simply the radial dis¬ 
tribution of matter, which can be directly obtained from 
the NFW profile. The function for the satellite-mass 
pairs, Fg^(x), can be derived by convolving satellite and 
dark matter profiles. Under our assumption that the 
satellite distribution inside a halo also follows an NFW 
profile, Fgj^(x) is the convolution of two NFW profiles 
with different concentration parameters, which has an 
analytic expression given in Appendix A. More general 
discussions of galaxy-galaxy leasing in the halo occupa¬ 
tion framework are given by Mandelbaum et al. (2005), 
who focus on halo virial masses and satellite galaxy frac¬ 
tions, and Yoo et al. (2006), who focus on the derivable 
constraints on flm and erg. 


3.4. The Galaxy Bispectrum 

At the final step of our analysis, we consider the large- 
scale behavior of the galaxy bispectrum, the Fourier 
transform of the three-point correlation function. Al¬ 
though this is a measure of real-space clustering, like 
those in §3.1, we treat it separately because its informa¬ 
tion content is rather different and because it is more 
difficult to measure observationally. The specific observ¬ 
able that we consider is the reduced bispectrum 


(3(ki, k 2 , ks) — 


i?(ki,k2,k3 


P{ki)P{k2) + P{k2)P{ko) + P{k3)P{ki) ’ 

(27) 

where i3(ki,k 2 ,kg) is the bispectrum, P{k) is the power 
spectrum, and the three wavevectors form a triangle 
(ki -|- k 2 -l- kg = 0). In the weakly nonlinear regime, the 
local relation between the galaxy density contrast 6g and 
the mass density contrast 6 can be written as a Taylor 
expansion Sg = ^bg^nd^/n\. The reduced galaxy bis¬ 
pectrum Qg is related to the reduced matter bispectrum 
Qm through 


(5g(ki,k2,k3) = 


Qm(kl,k2,k3) bg^2 
dg,l &g,l ’ 


(28) 


which shows contributions from both gravitational clus¬ 
tering and galaxy bias (Fry 1994; see also Fry & 


Gaztahaga 1993; Juszkiewicz et al. 1995). On large 
scales, 5gp and bg ^2 are simply the galaxy number- 
weighted halo bias factors b^p and bh, 2 , respectively, 
where bh,i and bh ,2 are predicted by nonlinear pertur¬ 
bation theory (Scoccimarro et al. 2001). The reduced 
matter bispectrum is calculated using the linear matter 
power spectrum P{k) and the second-order perturbative 
bispectrum (Fry 1984; Scoccimarro et al. 2001), 

BPT = 2 A 2 (ki,k 2 )P(fci)P(fc 2 ) +cyc., (29) 

where F 2 (ki,k 2 ) = j + i cosdi 2 (fci/A :2 + k 2 /k\) F 
|cos^ 0 i 2 , with ki • k 2 = kik2COs9i2- The weak de¬ 
pendence on Vim (Kamionkowski & Buchalter 1999) is 
neglected here. The key diagnostic power of Qg on large 
scales lies in the distinctive triangle-shape dependence of 
Qm predicted by gravitational perturbation theory. Pos¬ 
itive linear bias ( 6 gp > 1 ) suppresses Qg, and non-linear 
bias ( 6 g _2 > 0 ) can restore the amplitude but not the 
shape dependence (Fry 1994; Verde et al. 2002). 

4. BREAKING THE DEGENERACY BETWEEN BIAS AND 
COSMOLOGY 

Motivated by CMB anisotropy measurements (e.g., 
Netterfield et al. 2002; Pryke et al. 2002; Spergel et al. 
2003), the abundance of galaxy clusters (e.g., White, 
Efstathiou, & Frenk 1993; Eke, Cole, & Frenk 1996), 
high-redshift supernova observations (e.g., Riess et al. 
1998; Perlmutter et al. 1999; Riess et al. 2004), and 
measurements of the 3D galaxy power spectrum and 
the Lya forest power spectrum (e.g., Croft et al. 2002; 
Percival et al. 2002; Tegmark et al. 2004b; Seljak et al. 
2005b), we assume a central cosmology with Vim = 0.3, 
(jg = 0.9, Us = 1.0, and F = 0.2. We model a galaxy 
population of mean space density hg = O.Ol/i^Mpc”^, 
similar to that of galaxies with r-band absolute mag¬ 
nitude brighter than —19.5 -I- 51og/i in the SDSS. We 
therefore use Zehavi et al.’s (2005) results for Mr < 
—19.5-I-5 logh galaxies as a guide for our central model’s 
HOD parameters. We assume a satellite occupation 
S{M) = M/(1.4 X lO^^h“^M 0 ) in equation (10), i.e., 
the number of satellites is proportional to the halo mass 
well above the cutoff Mmin- We set the central galaxy 
cutoff parameter triogM to 0 . 2 , consistent with predic¬ 
tions from galaxy formation models (Zheng et al. 2005). 
With these choices, the cutoff mass required to match 
hg = O.Ol/i^Mpc”^ is Mmin = 5.55 x lO^h-iMg. The 
probability distribution of satellite galaxies is taken to be 
a Poisson distribution, i.e., w = 1 in equation (11). We 
assume that satellite galaxies have the same radial pro¬ 
file and velocity dispersion as dark matter within halos, 
i.e., A log Co = 0 and = 1 . 

For the central cosmology and the central HOD, we use 
the analytic formulae presented in § 3 to calculate a set 
of observables, which we thereafter treat as observational 
measurements. These observables are the rms galaxy 
overdensity (Jg(r) on scales of 8 h“^Mpc and 15/i“^Mpc; 
the group multiplicity function ngrp(> N) for N = 5, 
10, 20, and 40; the galaxy two-point correlation function 
^gg(r) at r = 0.1, 0.3, 0.5, 1, 2, 5, and 10 /i“^Mpc; the (3 
parameter; the galaxy pairwise velocity dispersion (T„(r) 
at r = 0.25, 0.5, 3, 5, and 10 h“^Mpc; the average virial 
mass of galaxy groups (Mvir(A^)) for N = 10, 20, 40, and 
80; the galaxy-mass cross-correlation function Vlmi^m{r) 
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at r = 0.25 and 0.5 h ^Mpc; and values of the reduced 
galaxy bispectrum on large scales, for wavevector trian¬ 
gles having two sides fixed at = 2fci = 0.05/iMpc“^ 
and subtending angles 0 = 0, 7r/4, 7r/2, 37r/4, and tt. We 
assume that each of the 30 “measurements” has a Ict frac¬ 
tional uncertainty of 10%, and that the measurement er¬ 
rors are uncorrelated. On its own, the 10% error assump¬ 
tion is probably pessimistic; even the largest scale observ¬ 
ables can probably be measured at least this well with 
the full SDSS, and some of the smaller scale quantities 
have already been measured more precisely (we investi¬ 
gate the effect of reducing errors for selected observables 
in §4.5). However, the assumption that the measure¬ 
ment errors are uncorrelated is optimistic, even though 
we space the observables for any given statistic (e.g., mul¬ 
tiple scales of the correlation function) fairly widely. Un¬ 
less there are strong positive correlations among many 
observables, the constraints on the HOD and cosmology 
should not degrade much with respect to our forecasts, 
but a full application to observations must include the 
covariances among observables, calculated from the data 
or from mock catalogs. 

Each time we change cosmological parameters, and 
thus the statistics of the halo population, we search for 
HOD parameters that yield the best possible match to 
the observables of the central model. Specifically, we 
minimize 


AX^=E 


{Ff - 


(30) 


where Fi represents the tth observable with an observa¬ 
tional error aFi^ the superscript “p” indicates the pre¬ 
dicted value for a given cosmology and HOD model, and 
the superscript “o” indicates the observed value calcu¬ 
lated for the central cosmological model and HOD pa¬ 
rameters. To isolate the information content of differ¬ 
ent clustering statistics, we start by considering only the 
spatial clustering observables described in §3.1, then add 
new observables one at a time. We use a Gauss-Newton 
scheme to carry out the minimization (see Appendix B 
for details). We characterize the acceptability of models 
in terms of Ay^, appropriate for “flat” priors that treat 
all parameter values equally. 


4.1. Constraining HOD Parameters for a Fixed 
Cosmology 

Before turning to cosmological parameter constraints, 
we first ask how well the HOD can be determined if we 
assume that the cosmological parameters are known per¬ 
fectly from independent data. Observational constraints 
on the HOD for an assumed cosmology have already been 
derived from the 2dFGRS, the SDSS, and other data 
sets, but these studies have all used restricted forms of 
the HOD and limited subsets of the galaxy clustering 
observables, such as the correlation function or the mul¬ 
tiplicity function (e.g., Jing & Borner 1998; Peacock & 
Smith 2000; Marinoni & Hudson 2002; Kochanek et al. 
2003; Lin et al. 2004; van den Bosch et al. 2003b; Maglioc- 
chetti & Porciani 2003; Zehavi et al. 2004, 2005; Collister 
& Lahav 2005). Here we adopt our much more flexi¬ 
ble HOD parameterization and examine the constraints 
that could be obtained with 10% measurements of all of 
the observables discussed above. Like Yan, Madgwick, & 
White (2003) and van den Bosch et al. (2005), we explore 


the HOD parameter space using a Monte Carlo Markov 
Chain (MCMC) technique (e.g., Gilks et al. 1996). We 
adopt the central cosmological model throughout, and we 
start the chain from the central HOD parameters, which 
by definition yield Ay^ = 0. At any point of the chain, 
we generate a new set of HOD parameters by taking a 
random walk in the HOD parameter space with the step 
size for each parameter drawn from a Gaussian distri¬ 
bution. We accept the new HOD model with a prob¬ 
ability of 1 if xLw < Xoid and exp[-(y2^^ - xlid)/‘^] if 
xLw > Xoid> where y^^^ and xlid are values of y^ for the 
new model and for the previous model, respectively. Flat 
priors are adopted for HOD parameters S{M), CTiogM 
and Lo in logarithmic space and A log cq and in linear 
space. 

The constraints on {N{M)) and other HOD parame¬ 
ters are illustrated in Figure 2. Figure 2a plots the en¬ 
velope of {N{M)) curves determined from HOD models 
that have Ay^ < 1, thus showing the l-cr uncertainty in 
{N{M)). With our adopted set of observables, the shape 
of {N{M)) near the cutoff mass is poorly constrained. 
The mean galaxy density imposes a strong constraint 
on the cutoff scale, but either hard or soft cutoffs can 
produce the same fig. The halo bias factor is nearly 
independent of mass in the low-mass regime (see, e.g., 
Jing 1998; Sheth, Mo, & Tormen 2001; Seljak & Warren 
2004), and these halos contain only a single galaxy, so 
variations in the cutoff profile do not affect either the 
two-halo or one-halo terms of ^gg(r’). Observables that 
we have not considered, like the void probability function 
or the Tully-Fisher relation, might provide stronger con¬ 
straints on the location and form of the {N{M)) cutoff. 

At higher halo masses, our observables provide much 
greater leverage, and {N{M)) is tightly constrained over 
the range (N) ~ 1.2 to (TV) ~ 100. The constraints 
loosen at still higher halo masses, since these halos are 
too rare to have much impact on the correlation function, 
and we only consider group virial masses up to A = 80 
and the multiplicity function up to A = 40. Figures 2b - 
2d show the marginalized likelihoods of the parameters 
to, A log Co, and with dotted vertical lines marking 
the Icr range. Under our assumption of 10% observa¬ 
tional errors, the distribution width parameter u can be 
constrained at the level of 10 — 15% (l-cr). For the dif¬ 
ference between galaxy and dark matter concentrations, 
the l-cr constraint is 0.10 — 0.15 dex, although with more 
small-scale points in fggir) this quantity could be pinned 
down more tightly. The 1-a uncertainty in the velocity 
bias is only a few percent, largely because the spa¬ 
tial clustering observables tightly constrain the true halo 
mass at fixed A, and the average virial masses then de¬ 
termine the value of 

Figure 3 shows correlations between some pairs of HOD 
parameters, inferred from the MCMC run. If galaxies 
are more centrally concentrated than dark matter inside 
halos, then the preferred width of P{N\M) narrows to 
suppress one-halo pairs in halos with small virial radii 
(Fig. 3a). A wider P(A|M) correlates, weakly, with a 
lower velocity bias (a„ < 1 for w > 1; Fig. 35), since 
more small-scale pairs involve satellite galaxies, which 
make larger contributions to the pairwise velocity dis¬ 
persion. There is a similar but stronger correlation with 
galaxy concentration (Fig. 3c), since positive Alogco al- 
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Fig. 2.— Uncertainties in the determination of the HOD with a known cosmology. Panel (a) shows the mean occupation functions of 
central galaxies {dashed line), satellite galaxies {dotted lines), and all galaxies {solid lines). Thick lines are for the input (central) HOD. 
Thin lines are envelopes of the mean occupation functions determined from HOD models with < 1, representing the 1-cr range. The 
other three panels show the marginalized likelihoods for the distribution width parameter lu, the halo concentration parameter Alogco, 
and the velocity bias respectively. Vertical dotted lines in each of these three panels mark the central 68.3% of the distribution (i.e., 
the l-(7 range). 


lows more small-scale pairs to come from massive halos 
with large velocity dispersions. Finally, if the mean oc¬ 
cupation at high halo masses increases (higher S '15 in 
Fig. 3d), then the mean halo mass at a given high mul¬ 
tiplicity decreases, so greater than unity is favored to 
keep the apparent virial mass (Mvir(iV)) near its original 
value. 

4.2. Changing ftm with as, Ug, and F Fixed 

Before considering general constraints in the 
(Om, CTs:r) parameter space in §4.5, we exam¬ 
ine constraints along several axes in this space: a pure 
change in , a pure change in as , “cluster-normalized” 
changes that preserve the value of as^L^ at fixed Ug 
and F, and linked changes in as, and F (at fixed 
Ug) that preserve the amplitude and slope of the halo 
mass function at high masses. As discussed by Zheng 


et al. ( 2002 ), the halo population responds to these 
changes in relatively simple ways, so we gain physical 
insight into the origin of the more general cosmological 
constraints by studying the behavior along these axes. 
These restricted constraints are also relevant to the case 
where other parameters are determined by independent 
data. 

Figure 4 illustrates the variation of the halo mass func¬ 
tion dn/d\ogM (Fig. 4a) and the non-linear matter cor¬ 
relation function ^mm(?') (Fig. 46) along these four axes. 
Open circles show the central model. Dotted lines show 
the effect of increasing ilm from 0.3 to 0.5, with other 
parameters fixed. Dashed lines show the effect of in¬ 
creasing as from 0.9 to 1.0. Dot-dashed lines show a 
cluster-normalized change, with Dm increased to 0.5 and 
as decreased to 0.7 to keep asCt^ fixed. Solid lines show 
the same (Dm, erg) combination, but with a correspond- 
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Fig. 3.— Marginalized contours and correlations between selected pairs of HOD parameters for a known cosmology. The parameters 
plotted are the distribution width parameter a;, the halo concentration parameter Alogcp, the velocity bias olv, and the mean occupation 
number S 15 in halos of mass Mq. In each panel, the plus sign indicates the values of the central HOD model, the two dotted 

contours correspond to Ax^ = 1 and 4, and the two solid contours are for Ax^ = 2.30 and 6.17 (68.3% and 95.4% confidence levels for two 
parameters). 


ing change in T that keeps the slope of dn/d\ogM fixed 
at the cluster mass scale 5 x (see Zheng et 

al. 2002). In all cases we use the Sheth & Tormen (1999) 
mass function formula (our eq [4]) and the Smith et al. 
(2003) prescription for computing the matter correlation 
function. Note that peculiar velocities of halos on large 
scales are proportional to similar to the scaling 

of velocity dispersions of cluster mass halos (Zheng et al. 

2002 ). 

The constraints for a pure Dm change are easiest to un¬ 
derstand, so we examine this case first and in the greatest 
detail. As shown in Figure 4, a pure Dm change simply 
shifts the mass function in proportion to Dm- The effect 
on the matter correlation function is very small, aris¬ 
ing entirely from changes in the concentrations of halos. 
Once halo masses are scaled in proportion to Dm, the spa¬ 
tial clustering of halos at a given scaled mass is virtually 
unchanged, although the velocities change in proportion 


to D))^® (Zheng et al. 2002). 

Figure 5 shows for the best-fit HOD model at each 
value of Dm, as different observables are added in suc¬ 
cession. Figures 6 and 7 show the parameters of these 
best-fit HODs. If we include only the spatial cluster¬ 
ing observables <Jg(r), ngrp(> V), and ^gg(r), we obtain 
results virtually identical to those of the central model 
by simply shifting the mass scale of the mean occupa¬ 
tion function in proportion to Dm- For these observ¬ 
ables alone, we get Ay^ = 0 at all Dm {horizontal dotted 
line in Fig. 5). Figure 6e plots the best-fit {N{M)) as 
a function of the scaled halo mass M x (0.3/Dm)- Be¬ 
cause the constraints imposed by our observables on the 
low end of {N{M)) are weak, the cutoff profiles of the 
scaled mass functions vary, but above (N) ~ 1.2 they 
are nearly indistinguishable. Thus, the spatial clustering 
observables effectively fix {N{M)) in terms of scaled halo 
masses, and only small further adjustments can be ac- 
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Fig. 4.— Variations of (a) halo mass function and (fc) mass 
correlation function with cosmology. The fitting formula of Sheth 
& Tormen (1999) for the halo mass function and that of Smith 
et al. (2003) for the nonlinear mass power spectrum are used for 
the calculations. In each panel, we show the results for the central 
cosmological model (Qm = 0.3, as = 0.9, and F = 0.2; open cir¬ 
cles), a pure change in Om to 0.5 {dotted line), a pure change in 
as to 1.0 (dashed line), a cluster-normalized model with Qrn = 0.5 
and (T8 = 0.7 (dot-dashed line), and a cluster-normalized model in 
which we also change the shape parameter to F = 0.11 to preserve 
the halo mass function as closely as possible (solid line). 

commodated when other observables are introduced. As 
shown in Figure 7, the best-fit models have a; = 1 just 
like the central model HOD, and they have Alogcp « 0, 
although there is a slight reduction in halo concentration 
at low Qrn to compensate for the higher concentrations 
of the halos themselves. 

The Qm degeneracy breaks as soon as we include dy¬ 
namical measures that are sensitive to the absolute mass 
or velocity scale of halos. The first such measure that we 
add is /3 = Q^/bg. For low Qm, the cutoff of {N{M)) 



Fig. 5.— Values of Ax^ as a function of flm, relative to the true 
model with flm = 0.3, in a cosmological model sequence where only 
flm is varied. As each set of observables is added in the analysis 
(in the order indicated in the legend), the for each cosmological 
model is minimized by solving the best-fit HOD parameters. Dif¬ 
ferent line (point) types show a sequence of Q.rn constraints from 
including more and more complementary observables. The inset 
box shows the parts of the lines with Ax^ < 6. 

smooths out so that more massive galaxies come from 
lower mass halos with lower bias factors [dashed lines in 
Figs. 6a and 66). The P[N\M) distribution broadens 
slightly (w > 1 in Fig. 7a) to compensate for a drop in 
one-halo pairs, and concentrations drop (Alogcp < 0 in 
Fig. 7b) so that these pairs move to larger separations. 
However, within the constraints imposed by the multi¬ 
plicity function, the correlation function, and the galaxy 
number density, there is virtually no room to change the 
galaxy bias factor in a way that compensates changes in 
Qm- For Qm in the range 0.2 — 0.4, the Ay^ values in 
Figure 5 are almost equal to those expected for constant 
bg, given our adopted 10% uncertainty in (3 (and hence 

m / 

When we add galaxy pairwise velocity dispersions as 
observables, the velocity bias parameter adjusts to 
compensate for changes in Qm (Fig. 7c, long-dashed line). 
However, while velocity bias can reduce velocity disper¬ 
sions within halos, it cannot affect the pairwise velocities 
of the halos themselves, which scale as Q^. For high Qm, 
drops by more than the factor (0.3/Hm)^^^ that would 
keep internal dispersions fixed, in an effort to compensate 
for the higher halo velocities. This change would drive 
the pairwise dispersion too low on small scales, where 
the one-halo term dominates, so halo concentrations in¬ 
crease, allowing more of the small separation pairs to 
come from massive halos. In response to these higher 
concentrations, the width parameter to drops, reducing 
the one-halo contribution to ^gg(r). Despite the freedom 
introduced by velocity bias, the combination of pairwise 
dispersions at small and large scales adds significant dis¬ 
criminatory power on the value of Qm- For Qm = 0.4, for 
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Fig. 6. — {a-d) Changes in {N{M)) as new sets of observables are added, for cosmological models differing only in Qm- The four 
panels are for flrn =0.1, 0.2, 0.4, and 0.5, respectively. For visual clarity, we omit the line showing the effect of adding the reduced galaxy 
bispectrum Qg. (e) Best-fit {N{M)) for different values of Qrn if only spatial clustering observables are considered, after scaling halo masses 
by (rim/O.S). (/) Same scaled mass result if other clustering observables are also included. Open circles in each of the six panels show 
{N(M)) for the central cosmological model. 


example, increases from ~ 3.5 {(3 alone) to ~ 5.6, 
with the two dynamical measurements producing most 
of the discrepancy. 

When we add the (Mvir(A)) observables, the velocity 
bias is forced close to the value = (0.3/Dm)^^^ that 
keeps these apparent virial masses fixed (Fig. 7c, dot- 
dashed line, nearly obscured by solid line). Since there is 
no longer room to improve the large-scale pairwise dis¬ 
persion by “overcompensating” with Uy, the concentra¬ 
tion and width parameters return to the values preferred 


by the spatial clustering observables (note the similarity 
of dotted and dot-dashed lines in Figs. 7a and 7b, and 
the difference between these lines and the long-dashed 
lines). Although ay can exactly compensate a pure flm 
shift for the (Mvir(A)) observables, the match to the 
pairwise dispersion becomes worse as a result, and 
increases slightly. 

The galaxy-galaxy lensing observable DmCgm(r) pro¬ 
vides a diagnostic of host halo masses that is unaffected 
by velocity bias. If {N{M)) stays fixed as a function 
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Fig. 7.— Other parameters {uj, A log cq, and ai,) of the best-fit 
HODs as a function of Qrn and the set of observables, for cosmo¬ 
logical models differing only in 

of scaled mass, then a pure Qm change simply scales 
in proportion to flm (eq. 25). The effect of 
higher flm can be partly compensated by reducing halo 
concentrations (Fig. 76, solid line), which dilutes the con¬ 
tribution from satellite galaxy and matter particle pairs. 
This decrease drives an increase in uj (Fig. 7a) to restore 
lost one-halo galaxy pairs at small separations (note that 
^gm itself is independent of w). However, altering con¬ 
centrations has limited power to compensate for a shift 
in the halo mass function, and Hm^gm('c) measurements 
add substantially to Ay;^. Measurements in the linear or 
near-linear regime, not considered here, would add still 
more discriminatory power. 

Finally, large-scale bispectrum measurements have lit¬ 
tle sensitivity to a pure shift (Fig. 5, long-dashed 
line). This statistic is mainly a diagnostic for the galaxy 
bias factor bg. Along this model sequence, the ampli¬ 
tude of dark matter clustering does not change, and the 
galaxy bias factor is therefore determined by the ob¬ 
served spatial clustering, independent of The bis¬ 
pectrum will become important when we consider se¬ 
quences with changing erg- 

Overall, these galaxy clustering observables have great 
power to constrain any pure shift in 0^- A model with 
Vljn ^ 0.3 has a likelihood exp(—Ax^/2) relative to the 
central model, so even a model with Ay^ = 4.6 is dis¬ 
favored by a factor of 10. Including all observables, we 
find Ax^ ~ 5.8 and 5.0 for 0^, = 0.25 and 0.35, respec- 



Fig. 8.— Values of as a function of erg, relative to the true 
model with erg = 0.9, in a cosmological model sequence where only 
(Tg is varied. The format is similar to Fig. 5. 

lively, implying that AHm = 0.05 can be firmly ruled 
out. These constraints weaken significantly if we omit 
the weak-lensing observables with their direct sensitivity 
to mass scales. Still, even with (3 and o'jj(r) as the only 
dynamical observables, the flm = 0.2 and 0.4 models are 
rejected with Ax^ > 5. 

4.3. Changing as with flm, ns, and F Fixed 

We now consider models that differ only in ug, the am¬ 
plitude of the linear matter power spectrum, with 
and the shape of P{k) held fixed. The dashed lines 
in Figure 4 show the effect of raising trg from 0.9 to 
1.0. The amplitude of the matter correlation function 
increases at all scales, and the space density of high-mass 
halos (M > \Q^'^h~^MQ) rises, while the space density 
of galaxy mass halos (M ~ 10^^ — is nearly 

unchanged. Figure 8 shows Ay^ for the best-fit HOD 
models along the changing erg sequence, and Figure 9 
shows the properties of these best-fit models. 

Figure 96 shows the response of {N{M)) to trg when 
the spatial clustering observables Ug, ^gg and rig^p are 
the only constraints. Higher as leads to lower (N{M)) 
at high masses, so that Ugrp remains fixed despite an in¬ 
creased halo abundance. This change reduces the galaxy 
bias factor bg , but on its own it is not enough to compen¬ 
sate for the increased matter clustering amplitude. The 
{N{M)) cutoff therefore spreads to lower masses, putting 
more isolated galaxies in low bias halos while retaining 
the galaxy number density Ug. The P{N\M) width pa¬ 
rameter rises to w > 1 to restore small-scale pairs lost 
when the number of halos in the {N{M)) ~ 1 — 2 range 
is reduced (Fig. 9c, dotted line). 

For spatial clustering observables alone, these HOD 
changes effectively compensate for changes in the halo 
population, and the Ay^ curve is nearly flat in the vicin¬ 
ity of the central model (Fig. 8, dotted line). However, 
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Fig. 9.— Best-fit HODs for cosmological models differing only in as- (a) (N{M)) curves for the as = 1.0 cosmological model as different 
observables are included, (b) Best-fit {N{M)) curves for a range of crs values when only spatial clustering observables are considered. The 
central model {open circles) has as = 0.9. (c, d, e) Other parameters {lj, Alogco, and a.v) of the best-fit HODs as a function of as for the 
same set of observables shown in (a). In (e), the dot-dashed line is obscured by the solid line. 


these observables pin down {N{M)) fairly precisely, leav¬ 
ing little room for adjustment when dynamically sensi¬ 
tive observables are added as constraints. The from 
P = is close to that expected for our 10 % error 

bar and the scaling bg oc l/us implied by fixed galaxy 
spatial clustering. The velocity bias parameter can 
partly compensate for changes in halo velocities and in¬ 
ternal dispersions, but it cannot fix the small- and large- 
scale values of cr„(r) simultaneously (see discussion in 
§4.2), so Ay^ rises when cr^ is added. Adding (Mvir(Af)) 
effectively pins down to the value needed to keep the 
apparent virial masses fixed, making the overall fit of 
av{r) worse and further increasing Ay^. The Dm^gm('c) 
observable is unaffected by velocity bias, and the small 
adjustments still allowed in {N{M)), Alogco, and w do 
little to counter the changes in halo masses. Finally, the 
galaxy bispectrum increases Ay^ through its direct sen¬ 
sitivity to bg. With all observables considered together, 
Ay^ exceeds 4.0 for Acrg = ±0.1, and it rises rapidly for 


larger changes. 

4.4. Cluster-normalized Models 

The observed abundance of massive galaxy clusters ef¬ 
fectively constrains the amplitude of the halo mass func¬ 
tion at a scale M ^ 5 x The halo abundance 

at this scale is an increasing function of flm and as, and 
the cluster normalization constraint is usually expressed 
in the form = const., with q ^ 0.5 (e.g., White, 

Efstathiou, & Frenk 1993; Eke, Cole, & Frenk 1996; Pier- 
paoli, Scott, & White 2001). Here, as in Zheng et al. 
( 2002 ), we consider a sequence of cluster-normalized cos¬ 
mological models that satisfy as = 0.9(f2m/0.3)“°-^, with 
the shape of P(k) held fixed. The dot-dashed lines in 
Figure 4 show the effect of simultaneously increasing flm 
to 0.5 and decreasing as to 0.70. As expected, the halo 
mass function stays the same at M ^ 5 x 
but relative to the fiducial model there are more low- 
mass halos and fewer halos with M > lO^®/i“^M 0 . The 
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Fig. 10.— Same as Fig. 5, but for a cluster-normalized model 
sequence where crgflJJ;® is held fixed as is varied. The corre¬ 
sponding values of erg are marked on the top axis. 

matter correlation function is lower at all scales because 
of the lower value of as- 

Figures 10 and 11 show the Ax^ and parameters of 
the best-fit HOD models along the cluster-normalized 
sequence, in the same format as Figures 8 and 9. The 
multiplicity function fixes the high-mass end of {N{M)) 
in a way that is nearly independent of flm, although there 
are slight variations in slope that reflect the changing 
slopes of the halo mass function (Fig. 11&). The higher 
abundance of low-mass halos for high requires an 
increase in the cutoff mass to keep the galaxy number 
density constant. The tightened mass range for single¬ 
galaxy halos leads to a larger one-halo contribution to 
^gg(r), which is compensated by decreasing the width 
parameter to and the concentration parameter Alogcp. 
Although the matter correlation amplitude is reduced for 
high flm and low erg, the bias of halos at fixed mass is 
higher, and these two effects conspire to keep the large- 
scale amplitude of the galaxy correlation function nearly 
constant over a substantial range of Qm- For spatial 
clustering observables alone, is nearly flat over the 
range 0.2 < flm < 0.5. 

Because cluster normalization keeps the mass and ve¬ 
locity scale of halos roughly fixed, the addition of dy¬ 
namically sensitive observables makes relatively little dif¬ 
ference to best-fit HOD parameters, and Ay^ increases 
more slowly than it does for pure flm or pure erg changes 
(note the smaller vertical scale of Figure 10 relative to 
Figures 5 and 8). The detailed changes in w, Alogco, 
and ay as new observables are added mostly reflect 
the changes in the cutoff of {N{M)) and corresponding 
changes in the fraction of single-occupancy halos. The 
Ay^ curves are asymmetric in flm because of the non¬ 
linear relation between and erg. For example, cluster 
normalization requires a 13% reduction in as to 0.78 at 
ftm = 0.4 but a 22% increase to erg = 1.10 at ftm = 0.2, 


and the still larger erg values required at lower are 
easily ruled out. The bispectrum constraint is especially 
useful in distinguishing cluster-normalized models be¬ 
cause it responds to bg independent of flm- Despite the 
compensating effects of as and flm in cluster-normalized 
models, the full set of observables yields Ay^ ~ 5.5 for 
flm = 0.4 and Ay^ ~ 13.1 for Qm = 0.2. 

As shown by Zheng et al. (2002), one can construct 
cluster-normalized model sequences that match the am¬ 
plitude and slope of the halo mass function at cluster 
scales by changing the power spectrum shape parameter 
in concert with and erg. The solid lines in Figure 4 
show a model with Qrn = 0.5, trg = 0.71, and F = 0.11. 
The halo mass function is nearly identical to the central 
model’s for M > 10^^’®/i“^Mq, although it is higher at 
low halo masses. The shape of the matter correlation 
function is significantly different because of the different 
shape of P{k). 

Figure 12 shows Ay^ as a function of flm along this 
model sequence. The required values of erg and F are 
shown in Figure 13 of Zheng et al. (2002), and the values 
of F are also marked on the top axis of our Figure 12. 
As speculated by Zheng et al. (2002), matching the slope 
of the halo mass function by changing the shape of P{k) 
does more harm than good: Ay^ values at a given flm are 
always higher than they are for the fixed P{k), cluster- 
normalized models shown in Figure 10. With all observ¬ 
ables included, models with flm = 0.2 and = 0.4 
have Ay^ ~ 17.6 and Ay^ ~ 10.2, respectively. 

4.5. More General Cases 

The four model sequences considered in Figures 5, 
8 , 10, and 12 trace four different curves through the 
(Dm, erg, F) space of cosmological parameters. We now 
examine the constraints and degeneracies in this pa¬ 
rameter space in more general terms. The computa¬ 
tional expense of finding the best-fit, 10-parameter HOD 
model for each cosmological model makes a comprehen¬ 
sive study of the full 3D space difficult. However, the fact 
that changing the power spectrum shape to match the 
halo mass function makes models easier to distinguish 
suggests that the “F” dimension is largely decoupled 
from the ~ erg plane. Furthermore, the galaxy power 
spectrum shape can be measured directly on large scales 
(e.g., Tegmark et al. 2004a; Cole et al. 2005; Tegmark 
et al. 2006), and we have not incorporated these direct 
constraints into our observables. 

To begin, therefore, we fix F = 0.2 and map out Ay^ 
contours in the — erg plane, finding the best-fit HOD 
parameters that minimize the values of Ay^ at each point 
on the cosmological parameter grid. We include all the 
observables discussed in §3 and incorporated in our ear¬ 
lier constraints. The solid contours in Figure 13a show 
Ay^ = 2.30 and Ay^ = 6.17, corresponding to 68.3% 
and 95.4% confidence intervals for two parameters. The 
most degenerate direction seems to be crgH^ = const, 
with q ^ 0.7 — 0.8, which deviates from the degener¬ 
acy axis of cluster normalization {q ~ 0.5) or large-scale 
redshift-space distortions {q ~ 0.6). This deviation re¬ 
flects the steeper scaling favored by the galaxy-galaxy 
leasing constraint, which has g ~ 1 even on these non¬ 
linear scales (Yoo et al. 2006), and by pairwise velocity 
dispersions with ay as an adjustable parameter. 

Dotted contours in this panel show Ay^ = 1 and 4, 
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Fig. 11.— Similar to Fig. 9, but for the cluster-normalized model sequence. 


from which one can read projected Itr and 2a constraints 
on the individual parameters as and Qm- The remain¬ 
ing panels show the ID likelihoods for these parameters 
and the combination crgD^, inferred by integrating the 
likelihood function exp(—Ax^/2) over each bin of the rel¬ 
evant quantity. Under our assumption of 10% fractional 
error in each observable, the l-cr uncertainty in the mat¬ 
ter density parameter flm is about 11%. The amplitude 
as of the matter power spectrum has a l-cr constraint 
that is better than 10%. The l-cr uncertainty in the best 
constrained quantity, crgD^ with q = 0.75, is at the 5% 
level. If we force q to 0.5, the precision is almost equally 
good. 

What constraints could we infer on T from the galaxy 
clustering measurements considered in our observable 
set? To address this question in approximate terms, we 
fix Um and erg to the central model values of 0.3 and 
0.9 and vary T alone, again finding best-fit HOD pa¬ 
rameters and minimum Ay^. We obtain Ay^ = 1 at 
T = 0.183, 0.219 and Ay^ = 4 at T = 0.167, 0.240. The 
10 % (Icr) constraint on T is similar to that found along 


our matched halo mass function sequence (Ay^ = 1 at 
T = 0.177, 0.223 and Ay^ = 4 at T = 0.161, 0.250; 
Figure 12), supporting our conjecture that the (flm,crs) 
constraints are mostly orthogonal to the P{k) shape con¬ 
straints. With our 30 observables and assumption of 
10 %, independent fractional errors, therefore, the full y^ 
surface approximately follows the contours of Figure 13 
with an orthogonal Gaussian of roughly 10% la width in 
the F direction. 

As discussed in the introduction, traditional ap¬ 
proaches to testing cosmological models with large-scale 
structure focus on the perturbative regime, where it is 
reasonable to adopt linear theory and linear bias for de¬ 
scribing the galaxy power spectrum and redshift-space 
distortions and second-order perturbation theory and a 
quadratic bias model for describing the galaxy bispec¬ 
trum. It is interesting to compare our HOD-based con¬ 
straints to those that could be derived from the same 
observables with the perturbative approach. Specifically, 
we use the Ug, /3, and Qg observables adopted previously, 
with the same 10% uncertainties. We employ a two- 
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Fig. 12.— Same as Fig. 5, but for a model sequence in which 
(Tg and F are varied as a function of to preserve the halo mass 
function as closely as possible. The corresponding values of the 
shape parameter F are marked at the top axis. 

parameter quadratic bias model, and we calculate the 
observables for a given cosmology and bias using ag = ba, 
(3 = /b, and Qg = Qm/h + where a and Qm 

are the rms fluctuations and reduced bispectrum for the 
dark matter. We examine an (O^, erg) grid with T = 0.2 
and find the values of b and &2 that minimize for 
each cosmological model. 

The top left panel of Figure 14 compares our HOD- 
based constraints (repeated from Figure 13) to those ob¬ 
tained from the perturbative analysis. The principal de¬ 
generacy axis for the latter constraints is = const., 

since /3 constrains /b and Og constrains bas- The 
bispectrum constraint on b breaks this degeneracy and 
produces closed contours. However, the constraints from 
the full HOD-based analysis of all observables are much 
tighter than the constraints from the smaller set of ob¬ 
servables that can be modeled by the perturbative ap¬ 
proach alone. Of course, the perturbative cosmology con¬ 
straints can be improved by including more observables, 
such as large-scale weak-lensing measurements or CMB 
anisotropies, but these observables will also improve the 
HOD-based constraints. Thus, despite the need for a 
more complex, multi-parameter model of galaxy bias, the 
ability of the HOD approach to model intermediate and 
small-scale clustering leads to substantial improvements 
in cosmological power. At the same time, an HOD anal¬ 
ysis of observed galaxy clustering yields valuable tests of 
theoretical models of galaxy formation. 

The obtainable constraints on cosmological parame¬ 
ters depend, of course, on the precision of the cluster¬ 
ing measurements. Our assumption of 10% fractional 
uncertainty on each observable is probably conservative, 
and current galaxy clustering data have already yielded 
higher precision for some of these observables, such as 
the two-point correlation function. To study the impact 


of improving measurement precision, we change the as¬ 
sumed error bars of different sets of observables by turns 
and recalculate contours. This investigation also 
highlights the relative cosmological sensitivity of the dif¬ 
ferent observables. Each panel of Figure 14 (except for 
the top left panel discussed previously) shows changes in 
the (HmjO-g) constraints caused by improving the mea¬ 
surement precision of one set of observables. The two 
dotted contours are the 1-cr and 2 -ct confidence levels 
(Ay^ = 2.30, 6.17) obtained by assuming 10% fractional 
errors for all sets of observables; they match the solid 
contours in Figure 13a. The solid contours are obtained 
by assuming 5% errors for the set of observables indicated 
in the upper right corner of the panel, keeping 10% for 
other sets. 

With the galaxy bias factor constrained by other ob¬ 
servables, a better measurement of galaxy density con¬ 
trasts fjg leads to a better constraint on erg, shrinking 
contours at both ends. The tighter ug constraint in turn 
produces a slightly tighter 0^ constraint, but there is 
essentially no improvement in the constraint (i.e., 

the width of the contours in the narrow direction is un¬ 
changed). Reducing the observational error bar on Qg 
produces a similar improvement for a different reason — 
the more precise measurements of Qg allow a better con¬ 
straint on the galaxy bias factor &g, so they help to break 
the degeneracy between hg and the amplitude of the dark 
matter power spectrum. 

Improving the pairwise velocity dispersion measure¬ 
ments produces modest improvements in both the crgH^ 
constraint and the individual constraints on erg and 0^- 
The former improvement comes mainly from the large- 
scale velocity dispersion, which depends on halo veloci¬ 
ties and thus scales roughly with The individ¬ 

ual parameter improvements come mainly from the in¬ 
terplay of small-scale with the apparent virial mass 
(Mvir(A)). Better measurements yield tighter con¬ 
straints on the velocity bias parameter, and the small- 
scale (T„ and (Mvir(A)) measurements then constrain the 
halo mass scale. Improving the (Mvir(A)) measurements 
on their own yields better erg and flm constraints but does 
not change the ers^^m constraint. 

Improving the Cgg(^) measurements leads to a signif¬ 
icantly better determination of erg for two reasons: it 
tightens the constraints on P{N\M) and thus on the 
large-scale galaxy bias factor, and it improves the mea¬ 
surement of the large-scale galaxy clustering amplitude 
itself. These two effects are analogous to those discussed 
earlier for Qg and cjg, respectively. Although ^gg(r) is 
not directly sensitive to 12^, the tighter erg constraint 
produces a tighter 12^ constraint when combined with 
dynamical measurements. Improving ngrp(> N) mea¬ 
surements produces surprisingly little improvement in 
the cosmological parameter constraints. If we instead 
degrade the fractional errors on ngrp(> N) to 50%, the 
constraints (shown by light solid contours) get noticeably 
but not dramatically worse. While it is always possible 
to reproduce ngrp(> N) exactly by choosing the appro¬ 
priate {N{M)), this choice pins down the mass associ¬ 
ated with a given multiplicity, and it largely pins down 
the galaxy bias factor. It appears, however, that the 
^gg(r) measurements already determine the high end of 
{N(M)) well enough for these purposes, so that moderate 
improvements in ngi.p(> N) do not add much power. If 
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Fig. 13.— Constraints on firn and erg (with F fixed) from galaxy clustering data, (a) Contour plot. The two dotted contours correspond to 
Ax^ = 1 and 4, and the two solid ones are for Ax^ = 2.30 and 6.17 (i.e., about the 68.3% and 95.4% confidence levels for two parameters). 
The two dashed lines represent {erg/0.9){Qm/0.3)'^ = 1 with q = 0.5 (thin) and q = 0.75 (thick). The other three panels show marginalized 
likelihoods on erg, Qm, and crgQfn, respectively. Vertical dotted lines in each of these three panels mark the central 68.3% of the distribution 
(i.e., the 1-cr range). In (d), thin and thick lines are for q = 0.5 and 0.75, respectively. 


we degrade the presumed measurement errors on Cgg(i’) 
to 25%, then the impact of improving the ngrp(> N) 
measurements is larger, but it is still not dramatic. The 
apparent redundancy of the ngrp(> N) and Cgg(i') con¬ 
straints makes them a valuable consistency check on the 
underlying assumptions of the HOD modeling. For ex¬ 
ample, the bias factor of single-galaxy halos affects ^gg(r) 
but not ngrp(> N), so an HOD model that makes the 
wrong assumption about this bias factor will lead to con¬ 
flicting estimates of {N{M)) at high masses. 

We have also examined the consequences of adopting 
priors on the HOD parameters, which could be motivated 
by theoretical or observational considerations. For exam¬ 
ple, numerical simulations give some idea of the plausible 
range of values for the velocity bias factor, and Tully- 
Fisher measurements can put limits on HOD parame¬ 
ters describing the cutoff profile of the mean occupation 
function. We characterize the priors by a l-cr deviation 


allowed from a central value and minimize the sum of 
the function and a quadratic term for the prior (see 
Appendix B). With a 1-a range of 1.00 ± 0.15 for the 
velocity bias factor or 1.00 ± 0.10 for the distribu¬ 
tion width parameter lu, we find that there is almost no 
improvement in cosmological constraints. The reason is 
that, even without these priors, the measurements have 
already excluded cosmological models that need extreme 
values of and to. For example, models that need 
as high as 1.3 have ^ 0.2 and as < 0.7 (see Fig¬ 

ures 7 and 9). Adopting a prior for the halo concentra¬ 
tion parameter Alogco with l-cr range of 0.0 ± 0.2 only 
slightly improves the cosmological constraints by exclud¬ 
ing some small regions of extremely high or low values of 
as- Priors on these HOD parameters must therefore be 
quite strong before they would improve the constraints 
obtainable from the observable set considered here. 

As we have shown in Figure 2, for a fixed cosmology. 
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Fig. 14.— Constraints on Vim and (js (with F fixed) from galaxy clustering data. Dotted contours in each panel are the same as the 
solid contours in Fig. 13, i.e., 68.3% and 95.4% confidence levels by assuming 10% fractional error in each of the 30 observables. Solid 
contours in the top-left panel are corresponding constraints using Ug, /3, and Qg and applying the perturbative regime model (see text). 
Solid contours in each of the other panels are constraints from HOD modeling by changing the fractional errors to be 5% for the set of 
observables indicated in the panel. In the panel for ngrp, we also show the case with the fractional errors in the multiplicity function 
changed to 50% (thin solid contours). 


the cutoff profile of the HOD cannot be well-constrained 
with the observables we adopt in this paper. This re¬ 
sult suggests that priors in the shape of the cutoff profile 
would not improve constraints on cosmological param¬ 
eters. To verify this expectation, we perform a test in 
which we fix the cutoff profile width criog m to the central 
model value. The resulting cosmological constraints are 
virtually identical to those in Figure 13, where CTiogM is 
not fixed. However, if the value of the absolute cutoff 
mass scale can also be pinned down, e.g. from the Tully- 
Fisher relation or other dynamical measurements in the 
single-galaxy regime, then cosmological constraints im¬ 
prove substantially. We have investigated the idealized 
case in which Mmin and aiog m are fixed to the values of 
the central model. We maintain the galaxy number den¬ 
sity by adjusting the amplitude of the satellite galaxy 
occupation function. The strong constraint on the halo 


mass scale from the value of Mmin produces a tight con¬ 
straint on rijn (0.300 ± 0.008), which in turn leads to a 
much tighter constraint on as (0.90 ± 0.04). 

Finally, we consider the constraints achievable with a 
brighter, lower density galaxy sample. Different sam¬ 
ples of galaxies have different HODs, but they should 
probe the same underlying cosmology. All the in¬ 
vestigations we have presented so far are based on 
the luminosity-threshold sample with a number den¬ 
sity hg = O.Ol/i^Mpc”^, roughly corresponding to the 
Mr < —19.5 + 51og/i galaxies in the SDSS. We increase 
the mass scale of our central model HOD to construct 
a galaxy sample with number density 10 times smaller 
{hg = O.OOl/i^Mpc”^), which approximates the SDSS 
Mr < —21 -I- 51og/i sample (Zehavi et al. 2005). With 
the same assumption of 10% fractional error on each of 
the 30 observables, we find that the cosmological con- 
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straints from these two galaxy samples are quite similar, 
with the degeneracy direction (TsO^ = const, slightly ro¬ 
tated to higher q for the lower density sample. Thus, 
measurements for these two largely independent galaxy 
samples should allow an important consistency check on 
cosmological conclusions and, if the results are combined, 
a ~ y/2 reduction in statistical errors. 

4.6. The “Influence Matrix” 

Perhaps the most important implication of Figure 14 
is that the cosmological constraints do not rely on one or 
two clustering statistics but instead emerge from the in¬ 
terlocking web of measurements. Figure 15 demonstrates 
this point in a different way. It shows the “influence ma¬ 
trix,” which we define in terms of the partial derivatives 
9 In Fi /9 In Oj evaluated at the central HOD and cosmo¬ 
logical model. Each predicted observable Fi is a function 
of cosmological and HOD parameters, and each Oj rep¬ 
resents one of these parameters, so the influence matrix 
encodes the response of the observables to the individual 
model parameters. When evaluating the partial deriva¬ 
tives, we adjust Mmin so that the galaxy number density 
stays fixed, but we keep the remaining HOD and cosmo¬ 
logical parameters (other than aj) fixed. The influence 
matrix is closely related to the Fisher matrix, 

daidaj \dai) \daj) ’ 

where F is the vector of the predicted observables and 
C is the covariance matrix of the errors. However, while 
the Fisher matrix depends strongly on the assumed er¬ 
ror covariance matrix C, the influence matrix is indepen¬ 
dent of this assumption, and it reveals the dependence of 
observables on model parameters in a more transparent 
manner. 

Figure 15 shows the influence matrix calculated at the 
central model. The area of each circle is proportional 
to the element of the matrix from the corresponding ob¬ 
servable {vertical axis) and parameter {horizontal axis). 
The HOD parameter Mmin, determined by matching the 
galaxy number density, is presented here as an observ¬ 
able, since its value could in principle be probed by 
dynamical or weak-lensing measurements in the single¬ 
galaxy regime. If two columns have similar elements, 
then the corresponding parameters are largely degener¬ 
ate, since an increase in one can be compensated by a 
decrease in the other. 

The first three columns of Figure 15 show how cluster¬ 
ing observables vary with cosmological parameters. Most 
elements of cHn F /d In and cHn F/cHn cts are positive. 
An increase in 12^ shifts the halo mass function to higher 
mass scales (Fig. 4), and with the fixed satellite HOD, 
this shift increases the group multiplicity function Ugrp. 
To keep a fixed galaxy number density, Mmin increases 
by a larger factor than Dm, and the shift of galaxies to 
halos with larger M/M* leads to stronger galaxy bias for 
spatial clustering statistics Ug and ^gg. An increase in 
(78 boosts the matter clustering directly. Raising either 
parameter increases halo masses and velocities, and thus 
the values of dynamical observables. The similar form 
of these vectors implies significant degeneracy between 
these two parameters, with an overall trend that roughly 
follows 91nF/cHnDm = 0.759InF/ cHucts, but the de¬ 
generacy is not complete. The power spectrum shape 


parameter F is largely decoupled from 0^ and cts- The 
much stronger sensitivity of Mmin to Dm than to any 
other cosmological or HOD parameter explains why a 
tight observational constraint on Mmin can produce such 
a tight constraint on 0^ (§4.5). The sensitivity of spa¬ 
tial clustering statistics to Ilm, which contrasts with the 
flat Ayfl shown by the dotted line in Figure 5, arises 
because we here keep all mass scales in the HOD fixed 
except for Mmin- If we instead scale all masses oc 0^, 
then the derivatives of ^gg, ngrp, and erg with respect to 
Ilm essentially vanish, as expected, while the derivatives 
of (Mvir(Af)) become one instead of zero. 

The rightmost nine columns in Figure 15 show how 
clustering observables vary with the HOD parameters, 
with cosmological parameters held fixed. A degener¬ 
acy between a cosmological parameter and HOD pa¬ 
rameters exists to the extent that the “influence vector” 
9 In F/9 In ac can be approximated by a linear combina¬ 
tion of the influence vectors 91nF/91na_f/^i, where ac is 
the cosmological parameter and an.i are the HOD pa¬ 
rameters. As expected from our earlier results, a sub¬ 
stantial degeneracy exists if we restrict F to the spatial 
clustering observables ^gg, Ugrp, and Cg. For example, an 
increase in as can be largely compensated by decreases 
in the parameters 5'i4 and Fis that describe the high end 
of {N{M)), and changes in uj and Alogco can take up 
some amount of residual difference. However, with the 
full set of observables, the structures of 91nF/91nDm, 
91nF/91n(T8) and 91nF/91nF are quite different from 
any of the influence vectors of individual HOD parame¬ 
ters. 

The visual structure of Figure 15 suggests that there 
is no linear combination of the HOD influence vectors 
that accurately approximates the cosmological parame¬ 
ter influence vectors. We verify this suggestion by solv¬ 
ing for the combination coefficients that minimize the 
squared difference between the influence vector and 
the corresponding linear combination of the HOD influ¬ 
ence vectors. (We do not include Mmin as an observable 
here, since it probably cannot be measured as precisely 
as the other quantities.) Figure 16a shows the HOD in¬ 
fluence vectors multiplied by these linear combination 
coefficients. The 5'i4 and Fis parameters dominate the 
best-fit linear combination, with compensating for the 
increase in (Mvir(Af)) that would otherwise occur. Fig¬ 
ure 16c shows that this best-fit linear combination {open 
squares) does not fully recover the flm influence vector 
{filled circles), with pairwise dispersions and (especially) 
galaxy-mass correlations providing the largest discrep¬ 
ancy. If as and F are allowed to vary in addition to the 
HOD parameters, then the best-fit linear combination 
{open cireles in Fig. 16c) is a better approximation to 
the flm influence vector, with clear improvement in ay 
and HmCgm- Figure 16& shows that the contribution from 
the as influence vector is much larger than that from any 
of the HOD influence vectors. Nonetheless, the changes 
in 5 'i4 and Fis significantly improve the agreement with 
the multiplicity function and small-scale correlation func¬ 
tion, relative to a pure as change with no freedom in the 
HOD. As argued earlier, the shape parameter F is largely 
decoupled from the strong cosmological degeneracy be¬ 
tween rim and as- 

Figure 15 also shows that clustering observables are in¬ 
sensitive to some HOD parameters, most notably An and 
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Fig. 15.— Logarithmic derivatives of the observables with respect to the cosmological and HOD parameters, starting from the central 
model. The area of each point is proportional to d\nF/d\na, where F is the observable listed on the vertical axis and a is the parameter 
listed on the horizontal axis. Open circles indicate positive derivatives and filled circles negative derivatives, and the scale is shown in the 
inset box. For observables with multiple scales, the smallest scale is shown at the bottom of the corresponding set of points and the largest 
scale at the top, and color coding to match circles to observables. In our modeling, the HOD parameter Mmin is determined by matching 
the galaxy number density once the other HOD parameters are set. We show it here as an “observable,” whose value could in principle be 
probed by dynamical or weak-lensing measurements. 


o'logMj which describe the low-mass end of the mean oc¬ 
cupation function. We reached similar conclusions from 
our analysis of {N{M)) constraints in §4.1 and the effect 
of a tJiogM prior in §4.5. The insensitivity of clustering 
observables to these parameters means that it is difficult 
to constrain them with clustering data. However, this 
insensitivity also means that poor knowledge of these pa¬ 
rameters does not introduce uncertainty in the cosmolog¬ 
ical conclusions. We have also investigated the influence 
matrix for a simpler HOD parameterization in which the 
mean satellite occupation function is a truncated power- 
law (Kravtsov et al. 2004; Zheng et al. 2005). The in¬ 
fluence vector of the power-law normalization is similar 
to that of Si 4 in Figure 15, while the influence vector of 
the power-law slope is similar in overall form but with 
a different mass and length-scale dependence, reflecting 
its differential impact on high- and low-mass halos. The 


influence of the low-mass cutoff is, again, weak. 

5. ENVIRONMENTAL VARIATION OF THE HOD 

The halo occupation function P{N\M) is, by defini¬ 
tion, a distribution weighted by halo number, the prob¬ 
ability that a randomly selected halo of mass M con¬ 
tains N galaxies. Our clustering calculations implicitly 
assume that this distribution and other parameters of the 
HOD like A log cq and are independent of the larger 
scale environment. A systematic dependence of a halo’s 
galaxy occupation on the overdensity 6 of its surround¬ 
ings would alter the clustering predicted for a given halo 
population and (number-averaged) P{N\M). Ignoring 
this dependence could therefore lead one to infer an in¬ 
correct HOD or incorrect cosmological parameters when 
fitting observed galaxy clustering. 

The assumption of an environment-independent HOD 
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Fig. 16. — Approximation of the flrn influence vector by a linear combination of other influence vectors. Panels (a) and (6) are in a format 
similar to Fig. 15, but each influence vector (each column) is multiplied by the corresponding best-fit linear combination coefficient (except 
for the flm inffuence vector in the first column). Panel (a) shows the case that keeps erg and F fixed and approximates the Qm inffuence 
vector by the HOD inffuence vectors only, and panel (6) shows the case that the Q,rn inffuence vector is approximated by the combination 
of the (Tg, F, and HOD inffuence vectors. Panel (c) shows the Q,rn inffuence vector {filled circles) and the best-fit linear combinations 
excluding the erg and F inffuence vectors {open squares) and including them {open circles). The open points track but do not completely 
overly the filled points, which indicates partial but not complete degeneracy between flrn and the other model parameters with respect to 
these observables. 


is rooted partly in the Bond et al. (1991) excursion set 
derivation of the extended Press-Schechter formalism. 
This derivation predicts that the statistical features of a 
halo’s assembly history depend only on its present mass, 
a point emphasized by White (1996). In agreement with 
this prediction, Lemson & Kauffmann (1999) and Sheth 
& Tormen (2004) found that fV-body halos of fixed mass 
in different environments have similar properties and for¬ 
mation histories, although Sheth & Tormen (2004) found 
a subtle trend for close pairs of halos to have higher for¬ 
mation redshifts. The resolution of the simulations used 
in these studies effectively limited them to halo masses 
M > IO^^Mq. Two recent studies (Gao et al. 2005; 
Barker et al. 2006), based on the higher resolution, Mil¬ 
lennium Run V-body simulation (Springel et al. 2005), 
show that the correlation between halo formation time 


and large-scale environment becomes much stronger for 
halo masses M < IO^^'^Mq. To the extent that galaxy 
properties depend on halo formation time and halo mass, 
this correlation will in turn produce an environmental 
dependence oi P{N\M). 

There are several reasons for thinking that the quan¬ 
titative impact of such a dependence is small, at least 
for galaxy samples defined by simple thresholds in lu¬ 
minosity or baryonic mass. Empirically, Blanton et al. 
(2006) show that the correlation of SDSS galaxy prop¬ 
erties with galaxy overdensity on a scale of 8/i“^Mpc 
can be entirely explained by the correlation with over¬ 
density on the lh~^ Mpe scale characteristic of large 
halos; at fixed local overdensity, they detect no resid¬ 
ual correlation with the larger scale environment. On 
the theoretical side, Berlind et al. (2003) show that the 
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HOD (more precisely, the mean occupation function) 
of galaxies above a baryonic mass threshold in Wein¬ 
berg et al.’s (2004) cosmological hydrodynamic simula¬ 
tion is independent of environment within the statistical 
uncertainties imposed by the (50/i“^Mpc)^ simulation 
volume. Yoo et al. (2006), analyzing the same simula¬ 
tion, show that explicitly eliminating any environmen¬ 
tal dependence by randomly shuffling the galaxy popu¬ 
lations among halos of similar mass changes the galaxy- 
galaxy and galaxy-matter correlation functions by < 5% 
on scales 0.1h~^ Mpc < r < 5h~^ Mpc, again within the 
statistical uncertainties. In a similar experiment using 
semi-analytic galaxy populations in the Millennium Run 
simulation, Croton et al. (2006) find that the correlation 
functions of luminosity-bin samples change by a few per¬ 
cent; because of the large simulation volume, these shifts 
are measured with high precision. 

Since we hope in the long term to achieve few per¬ 
cent precision on cosmological parameter measurements, 
the Millennium Run results imply that environmental 
dependence of the HOD must eventually be considered 
as part of the program outlined in this paper. Allowing 
arbitrary dependence of P{N\M) on large-scale overden¬ 
sity S would probably make the problem intractable and 
degenerate, but given the small magnitude of the antici¬ 
pated effects, a more restricted parameterization should 
suffice. For example, since the age dependence of halo 
clustering is found mainly in the mass range of individ¬ 
ual galaxy halos, it may well be enough to allow 
to vary with S, or to treat the effective bias parameter 
of low-mass halos as a parameter to be fitted instead of 
using the standard Wbody result. Since much of the sig¬ 
nal in most clustering statistics comes from multi-galaxy 
halos, addition of such parameters may not have much 
impact on the achievable constraints on cosmological pa¬ 
rameters. Semi-analytic mock catalogs will be valuable 
in guiding the choice of parameterizations and in testing 
whether fitting methods that adopt a particular parame¬ 
terization (or that ignore environmental variations alto¬ 
gether) yield incorrect estimates of HOD or cosmological 
parameters. 

For cosmological applications, it will be desirable to 
work with galaxy samples defined by properties that are 
sensitive to halo mass and relatively insensitive to forma¬ 
tion history, such as stellar mass or luminosity in a band 
dominated by old stellar populations. Color-selected or 
morphologically selected samples will be more suscepti¬ 
ble to environmental dependence of halo formation his¬ 
tory at fixed mass. Conversely, samples with high lumi¬ 
nosity or mass thresholds may be the least susceptible, 
since the correlation of formation time with environment 
is weak for high-mass halos. While the breakdown of the 
Bond et al. (1991) prediction is an annoyance for the cos¬ 
mological applications that are the focus of this paper, it 
could be an asset for efforts to understand the physical 
processes that determine galaxy morphology and spec¬ 
tral properties. For example, if fitting the clustering of a 
selected class of galaxies requires environmental depen¬ 
dence of the HOD, then the form and magnitude of that 
dependence will provide strong clues to the aspects of 
halo formation that determine whether galaxies are mem¬ 
bers of that class. Conversely, if the color and morphol¬ 
ogy dependence of observed clustering can be entirely 
explained without environmental variations of the HOD, 


it will imply that aspects of halo formation that correlate 
with environment do not play a strong role in determin¬ 
ing these properties of galaxies. 

6. SUMMARY AND DISCUSSION 

Our analysis shows that HOD modeling can substan¬ 
tially increase the cosmological power of galaxy cluster¬ 
ing measurements, by breaking degeneracies between the 
clustering of dark matter and the bias of galaxies with 
respect to mass. Changing the shape or amplitude of 
the matter power spectrum or the value of 0^ alters 
the mass function, spatial clustering, and velocity statis¬ 
tics of the dark halo population in well-understood ways 
(Zheng et al. 2002). Our experiments here, which verify 
the qualitative arguments of Berlind & Weinberg (2002) 
and Zheng et al. (2002), show that changes to the galaxy 
HOD cannot mask these changes in the underlying dark 
halo population. With our highly flexible parameteri¬ 
zation of the HOD, the set of observables considered 
here yields Ict uncertainties of ~ 10% in erg, 17^, and 
F and ~ 5% uncertainty in the combination crgD^ with 
q ~ 0.75. We expect these forecasts to be conservative, 
as we have not included observables for which we did not 
have ready analytic approximations, and our assumption 
of 10% measurement errors is pessimistic in at least some 
cases. 

The physical origin of these cosmological constraints 
is straightforward to understand for simple changes in 
Dm, CTg, or CTgDm®, as discussed in §§4.2-4.4. The general 
theme of these discussions is that, for a given cosmo¬ 
logical model, the spatial clustering of galaxies largely 
determines the number of galaxies in halos of a given 
spatial abundance. Dynamically sensitive statistics then 
reveal the halo mass scale, which depends on cg and 
Dm. We allow an arbitrary bias between the galaxy 
and dark matter velocity dispersions within halos, but 
this freedom does not eliminate the constraining power 
of dynamical observables because the space velocities of 
the halos themselves do not change. The main param¬ 
eter degeneracy is approximately because fixing 

this combination roughly fixes the halo velocity scale and 
the abundance of halos at the mass scale of rich galaxy 
groups. However, the changing shape of the mass func¬ 
tion, the differing sensitivities of different velocity mea¬ 
sures, and the different (CTg,Dm) dependence of galaxy- 
galaxy lensing all serve to break this degeneracy. Fig¬ 
ures 14 and 15 demonstrate that the cosmological con¬ 
straints emerge from the full web of clustering observ¬ 
ables and are not dominated by one or two statistics on 
their own. 

Constraints on the galaxy HOD will themselves pro¬ 
vide valuable tests of galaxy formation models. The 
cutoff regime of {N{M)) is difhcult to pin down with 
the clustering statistics considered here, but for a known 
cosmological model the relation between average satellite 
number and halo mass is well determined, and the satel¬ 
lite distribution width oj, concentration Cg, and velocity 
bias are measured to ~ 10%, 30%, and 3%, respec¬ 
tively (Fig. 2). All of these quantities depend in detail 
on the physics that governs the evolution of satellites 
in larger halos (see, e.g., Taylor & Babul 2004, 2005a,b; 
Zentner et al. 2005), while the relative mass of halos that 
host central and satellite galaxies depends on the effi¬ 
ciency with which halos feed baryonic mass to their cen- 
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tral objects (see, e.g, discussions by Berlind et al. 2003; 
Zheng et al. 2005). These galaxy formation constraints 
will be especially powerful when derived as a function of 
luminosity, stellar mass, or other observables. 

The two main assumptions built into our modeling are 
the central-satellite parameterization and environment 
independence of the HOD. The central-satellite distinc¬ 
tion appears well rooted in galaxy formation physics, and 
it allows us to represent the range of plausible galaxy 
HODs more completely with a moderate number of pa¬ 
rameters. However, we have confirmed with other tests 
that if we model galaxy bias with a flexible HOD param¬ 
eterization that does not impose a central-satellite dis¬ 
tinction, but instead introduces a characteristic mass for 
the narrow-to-wide transition of P{N\M) as in Scranton 
(2003), then we reach almost identical conclusions about 
the cosmological constraining power of the clustering ob¬ 
servables considered in this paper. Recent numerical re¬ 
sults imply an environmental dependence of halo forma¬ 
tion times that opens the door to environmental varia¬ 
tion of the HOD, especially in the single-galaxy regime. 
As discussed in §5, we expect the quantitative impact of 
such dependence to be small, but potentially significant 
at the high level of precision we ultimately hope to at¬ 
tain. Investigation of environmental dependence effects 
and methods of allowing for them in HOD modeling are 
a high priority for future work. 

The cosmological modeling approach advocated here is 
closely related to the CLF method introduced by Yang, 
Mo, & van den Bosch (2003) and van den Bosch, Mo, 
& Yang (2003a), who use clustering data and the global 
galaxy luminosity function to constrain the dependence 
of the luminosity function on halo mass. In principle, 
the CLF and HOD methods are equivalent — they are 
merely differential and integral forms of one another. 
One can derive the CLF from a series of HOD fits to 
galaxy samples with different luminosity thresholds (Ze- 
havi et al. 2005; Tinker et al. 2005). Conversely, one can 
integrate the CLF to infer {N{M)) for galaxies above a 
luminosity threshold (Yang, Mo, & van den Bosch 2003; 
van den Bosch et al. 2003b, 2005). The principal virtue 
of our HOD-based approach is that by focusing on a sin¬ 
gle, well-defined class of galaxies, we can parameterize 
the HOD in a way that seems likely to capture the pre¬ 
dictions of any reasonable galaxy formation model. This 
kind of comprehensive parameterization is more difficult 
to achieve for the full CLF, and most analyses to date 
have assumed, for example, that the CLF has a Schechter 
form in halos of fixed mass. Nonetheless, it is valuable 
to pursue both HOD and CLF approaches and test for 
consistency of conclusions. 

HOD modeling complements rather than replaces the 
perturbative approach based on large-scale measures 
that can be modeled with linear or quadratic bias. HOD 
modeling is more complex, but it can take advantage of 
high-precision clustering measurements on small and in¬ 
termediate scales. HOD modeling can also amplify the 
power of the perturbative approach, extending its reach 
further into the non-linear regime and checking its range 
of validity at a desired level of precision. For example. 
Tinker et al. (2006) show that an HOD-based approach 
to redshift-space distortions can improve recovery of the 
perturbative parameter /3 that controls large-scale flows. 
Yoo et al. (2006) show that the linear bias model pro¬ 


vides an accurate description of galaxy-galaxy leasing for 
r > 2/i“^Mpc, and they show how to accurately model 
this phenomenon on smaller scales (see also Guzik & Sel- 
jak 2002; Mandelbaum et al. 2005). J. Yoo et al. (2007, 
in preparation) show that the scale-dependent bias fac¬ 
tors derived by fitting the projected galaxy correlation 
function can extend recovery of the shape of the linear 
matter power spectrum into the mildly non-linear regime. 

CLF and HOD analyses of the 2dFGRS and SDSS red- 
shift surveys have already produced a number of inter¬ 
esting results, even though they have considered only a 
fraction of the potential clustering observables. These 
results confirm, at least qualitatively, many of the basic 
predictions of current galaxy formation models, includ¬ 
ing the general form of the mean occupation function, 
the dependence of this function on luminosity, the exis¬ 
tence of a minimum mass-to-light ratio in the halos of 
~ L* galaxies where galaxy formation is most efficient, 
the large gap between the minimum halo mass for cen¬ 
tral and satellite galaxies above a luminosity threshold, 
the sub-Poisson fluctuations of P{N\M) that are a con¬ 
sequence of this gap, and the strong preference of galax¬ 
ies with older stellar populations for higher mass halos 
(van den Bosch, Mo, & Yang 2003a; van den Bosch et 
al. 2003b; Magliocchetti & Porciani 2003; Tinker et al. 
2005; Yang et al. 2005; Zehavi et al. 2004, 2005; Gollister 
& Lahav 2005). In combination with GMB data, the cos¬ 
mological constraints from HOD modeling of the SDSS 
projected galaxy correlation function are almost as tight 
as those from the large-scale galaxy power spectrum, and 
the two analyses are consistent within their statistical un¬ 
certainties (Abazajian et al. 2005). For the most part, 
the cosmological inferences from GLF/HOD modeling of 
galaxy clustering agree with those from other methods, 
but matching the mass-to-light ratios of galaxy clusters 
simultaneously with other clustering data appears to re¬ 
quire values of 0^ and/or erg that are substantially lower 
than the commonly adopted values of 0.3 and 0.9 (van 
den Bosch, Mo, & Yang 2003a; Tinker et al. 2005). If 
this conclusion is correct, then the evidence for it should 
become much stronger as more clustering observables are 
brought into play and the SDSS data set itself moves to 
completion.® 

We have focused in this paper on the cosmological pa¬ 
rameter constraints that can be derived from galaxy clus¬ 
tering alone, using external data only to guide the choice 
of power spectrum shape and motivate the assumption of 
Gaussian initial conditions. As with perturbative analy¬ 
ses of large-scale structure, the long-term interest lies in 
combining these constraints with those from the CMB, 
the Lya forest. Type la supernovae, and other cosmo¬ 
logical observables. The complementary sensitivities of 
these observables lead to much tighter parameter con¬ 
straints. More importantly, conflicts among them could 
point the way to physics beyond the simplest versions 
of ACDM, such as evolving dark energy, a gravitational 
wave contribution to GMB anisotropy, departures from 
scale invariance in the primordial power spectrum, non- 

® Since the original submission of this paper, the analysis of 
the three-year WMAP data set (Spergel et al. 2006) has provided 
strong support for this shift in cosmological parameter values, in 
excellent agreement with the results of the HOD and CLF model¬ 
ing and the related “empirical model” method of Vale & Ostriker 
(2006). 
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zero space curvature, cosmologically significant neutrino 
masses, and so forth. By sharpening the constraints from 
large-scale structure in the new generation of galaxy red- 
shift surveys, HOD modeling can play a critical role in 
efforts to test the standard cosmological model and, per¬ 
haps, discover its breaking points. 
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APPENDIX 

PAIR DISTRIBUTIONS OF TWO POPULATIONS FOLLOWING NFW PROFILES WITH DIFFERENT CONCENTRATIONS 


For the calculation of the one-halo term of the galaxy-galaxy or galaxy-mass two-point correlation function, we need 
to know the distribution of galaxy pairs or galaxy-mass particle pairs as a function of separation in each halo. The pair 
distribution is basically the convolution of the density profiles of two populations. For spherically symmetric density 
distributions pi and p 2 , the fraction of pairs with separation in the range of r to r -|- dr is 


dF f 

—droc / pi(xi)d^xi p 2 {x 2 +r)d^r 

oc J 4:TrxiPi{xi)dxi j p 2 {x 2 )‘ 2 ,'xr^ sinadadr 

^+oo 

(xr^dr J dxixipi{xi) J dl3p2{x2), 


(Al) 


where F{r) is the fraction of pairs with separation less than r, = xf -I- r^ — 2xir/3, (3 = — cos a, and a is the angle 
between xi and r. If pi and p 2 follow the same NFW profile truncated at the virial radius, the analytic expression 
of dF/dr can be found in Sheth & Diaferio (2001) (see their eq.[l] and eq.[A25]; note that A(r) in their equations 
is proportional to r~^dF/dr). In this paper we assume that, like the dark matter, the distribution of galaxies inside 
halos also follows the NFW profile. Therefore, we use the formula in Sheth & Diaferio (2001) when evaluating the 
one-halo term of galaxy-galaxy correlations. However, we allow galaxies to follow a different NFW profile than the 
dark matter in halos of a given mass, which is our way to parameterize the spatial bias inside the halo in this paper. 
So, for the purpose of calculating the one-halo term of the galaxy-mass two-point correlation function, we need the 
pair distribution function corresponding to two NFW profiles with different concentration parameters. We evaluate 
the convolution in equation Al for two NFW profiles with concentration parameters ci and C 2 , both truncated at the 
virial radius i?vir of the halo. After tedious algebra, we derive the following analytic expression for the (differential) 
pair distribution function dF{x; ci, C 2 )/dx, where x = r/ (2i?vir)- 
For 0 < X < 0.5 


where 


dF{x;ci,C2) 

dx 

/i(x;ci,C 2 ) = 

/2(x;ci,C2) = 


= A(ci,C2)[/i(x;Ci,C2) -I- f2ix-,Ci,C2) + f3{x]Ci,C2) + fiix] Cl, C2)]s, 
1 


(C2 +Ci+ C1C2S)^ 


ln[(l -I- Cis)(l -I- C2s)] 


CiS 


C2(C2 +Ci+ C1C2S)(1 J- Cis) ' 


1 


/3(x;ci,C2) = 


(C2 - Ci+ CiC2s)^ 
1 


In 


(A2) 

(A3) 
(A4) 
(A5) 

C2(l + Cl)(l -I- C2)(l -I- Cl - Cls) ’ ^ 

and s = r/Ryir = 2x. For the purpose of stableness, one should compute /2 using the limit value —C 2 ^(l-l-ci)“^]/2 
as s approaches (ci — C 2 )/(ciC 2 ) and similarly /s = [cj"^(l + 02 )“^ — ]/2 as s approaches (c 2 — ci)/(ciC 2 ). 

dF(x;ci,C2) 


(C2 - Cl - CiC2s)2 

/4(x;ci,C2) = -- 


[ (1 -F Cis)(l -F C 2 - C 2 S) 


Cl(l - s) 

[ 

1 + Cl J 


C 2 (C 2 — Cl -F CiC 2 s)(l -F Cis)(l -F Cl) ’ 

- In 

(1 -F C 2 s)(l -F Cl - cis) 

1 Ci(l-s) 

2 

1 -F C2 

C2(C2 — Cl — CiC 2 s)(l -F Cl — Cis) 


For 0.5 < X < 1, 


dx 


= ^(ci,C2)/(x;ci,C2)s, 


where 


/(x;ci,C 2 ) = 


(C2 -I- Cl -I- CiC2s)2 


In 


(1 + Ci)(l J- C2) 


-F 


s- 2 


(1 — Cl + Cis)(l — C 2 + C 25 ) J (1 + Ci)(l + C2)(c2 + Cl + C 1 C 25 ) 


(A7) 

(AS) 
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The distribution is normalized so that fgdx dF/dx = 1. The normalization factor yl(ci,C 2 ) in equations (A2) and 
(A7) can be expressed as 

A(ci,C2) = \/Ah,(ci)Ah,(c2), (A9) 

where we have the following fitting formula for Ah<(c), 

Ah.(c) = Aoc^+“[ 1 + + Bosin[w(logc-()i)]}, (AlO) 

with Ao = 2.4575, a = -3.099, (3 = 0.617, ct = 1.651, = 4.706, Bq = 0.0336, w = 2.684, and 4) = 0.4079. The fitting 

formula has a fractional error less than 0.2% for 1 < c < 100. 

It can be shown that the pair distribution derived here (eqs.[A2] and [A7]) reduces to equation (A25) in Sheth & 
Diaferio (2001) for the case ci = C 2 . We also test the above formulae by generating random particle distributions 
following different NFW profiles and counting particle pairs directly. The Monte Carlo result agrees with the analytic 
formulae perfectly. 


AND MINIMIZATION 

We use the Gauss-Newton method to perform the minimization. Here we briefly describe the main procedure. 
Readers are referred to Gould (2003) for more details. 

Assume that we have n observables Fk (k = 1, 2,..., n) predicted by a model with m parameters {i = 1, 2,..., to). 
Observations give Fk^ohs- The is defined as (Einstein’s summation convention is used) 


y^ = [Tfc(a) — Fk^ohs]Bki[Fi{a) — Ei^obs], 


(Bl) 


where B^i is the inverse of the covariance matrix and has only diagonal components Bkk = if observational 

errors are uncorrelated. 

If the initial guess of the solution of a is a,, through linearizing E'fc(a) around a* and demanding that the first 
derivatives of y^ are zero, we obtain an equation for the correction to the initial guess (Gould 2003), 


bAa = d. 


where 


and 


bij 


dFk „ dFi 

dai ^ daj ^ 


1 d'^x^ 

2 daidaj ^ 


dt 


m 

dai 


Bki \Fi ,obs 


E;(a,)]. 


(B2) 

(B3) 

(B4) 


The direction of correction for a, is given by Aa. We set a-|-eAa as the new a,, where e is a small number, and refine 
the estimation of the solution of a. The above process is repeated until y^ reaches a minimum. It can be proved that 
the covariance matrix of a is just (Gould 2003). 

If some parameters have additional constraints, we introduce a cost function to realize the constraints and minimize 
the sum of y^ and the cost function. For example, if parameter oi cannot exceed 02 , the cost function can be defined 
as C = [(oi — a 2 )/cr]^ for oi > 02 and C = 0 for oi < 02 , where a is the tolerance. Another example is that a parameter 
has some priors. If oi is assigned a prior of Gaussian distribution with mean di and standard deviation cr, the cost 
function can be defined as C = [(oi — ai)/cr]^. Including a cost function leads to additional terms, —9C/9ai/2, in di 
(eq. [B4]). 
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